Internal Thermal Analysis of

Electric Machines
Improved Thermal Modelling of an
Electric Machine using Computational Fluid Dynamics
Master’s Thesis in Automotive Engineering

ANTHONY JAYANATH VIVEK

Department of Mechanics and Maritime Sciences

C HALMERS U NIVERSITY OF T ECHNOLOGY
Gothenburg, Sweden 2019
 Master’s thesis 2019:88

Internal Thermal Analysis of

Electric Machines

Improved Thermal Modelling of an

Electric Machine using Computational Fluid Dynamics

ANTHONY JAYANATH VIVEK

Department of Department of Mechanics and Maritime Sciences

Division of Vehicle Engineering and Autonomous Systems
Chalmers University of Technology
Gothenburg, Sweden 2019
Internal Thermal Analysis of Electric Machines
ANTHONY JAYANATH VIVEK

© ANTHONY JAYANATH VIVEK, 2019.

Supervisor & Examiner: Alexey Vdovin, PhD.,

Post Doctoral Researcher,
Department of Mechanics and Maritime Sciences,
Chalmers University of Technology

Master’s Thesis 2019:88

Department of Mechanics and Maritime Sciences
Division of Vehicle Engineering and Autonomous Systems
Chalmers University of Technology
SE-412 96 Gothenburg
Telephone +46 31 772 1000

Cover: Temperature distribution in the ERAD motor overlayed with a picture of

the ERAD motor used for simulations in this thesis[1].

Typeset in LATEX
Gothenburg, Sweden 2019

iv
Internal Thermal Analysis of Electric Machines
ANTHONY JAYANATH VIVEK
Department of Mechanics and Maritime Sciences
Chalmers University of Technology

Abstract
Electric machine thermal modelling using 3D Computational Fluid Dynamic (CFD)
simulations, is a vital process in the design and verification of Lumped Parameter
Thermal Networks (LPN). In this thesis, the aim is to continue the CFD Modelling
of a Permanent Magnet Synchronous Motor, and further improve the existing LPN
model. This thesis is a continuation of the Master’s thesis “Electric motor internal
heat convection modelling and analysis” conducted by Anil Kumar during 2018. The
CFD model is built using a bottom-up approach where each significant region of the
motor (Ex: Air-gap, End windings, etc) is modelled and validated by similar models
from journals, empirical formulation, literature or experiments. The complete CFD
model is then validated against rig measurements. With this assertion, the CFD
model is then used to run steady state simulations for various load cases and operat-
ing points to extract data. This data can be used as parameters in the LPN model,
or to calibrate the model itself. This way, the highly simplified 1-Dimensional LPN
model is verified by 2 or 3-Dimensional CFD simulations, which itself is validated
by literature and physical measurements. These 1D LPN, 2D models or 3D CFD
simulations can be applied in the future, for various different electric machines or
prototype models to analyze cooling strategies and performance.

Keywords: Computational Fluid Dynamics, Electromagnetic, Electric Machine,

PMSM, Rotor, Stator, End windings, Thermal, LPTN, Drive Cycle Analysis.

v
Abstrakt
Termisk modellering av elektriska maskiner med hjälp av tredimensionella numeriska
beräkningsmetoder inom CFD (Computational Fluid Dynamic) är en viktig aspekt
vid design och verifiering av termiska LPN-modeller (Lumped Parameter Networks).
Denna studie syftar till att numeriskt modellera en permanentmagnetiserad synkro-
nmotor och vidareutveckla en redan existerande LPN-modell. Detta examensarbete
är en fortsättning på Anil Kumars examensarbete ”Electric motor internal heat
convection modelling and analysis” från 2018. Den numeriska CFD-modellen är
uppbyggd genom en bottom-up-metod, där varje väsentlig del av motorn (exem-
pelvis luftgapet och lindningarna) modelleras och valideras med hjälp av enklare
modeller från vetenskapliga artiklar, empiriska formuleringar och litteratur. Den
fullständiga numeriska CFD-modellen valideras med hjälp av experimentella data.
Stationära simuleringar utförs med den numeriska modellen för olika lastfall och
driftspunkter. Data från dessa simuleringar används som parametrar för att kalibr-
era LPN-modellen. Genom användning av de tredimensionella simuleringarna valid-
eras den mycket förenklade, endimensionella LPN-modellen, vilka i sin tur valideras
av resultat i vetenskapliga artiklar och experimentella mätningar. Såväl den förbät-
trade CFD-modellen som LPN-modellen är lämpliga för att genomföra analysera
kylstrategier och kylprestanda för olika elektriska maskiner och prototyper.

Nyckelord: Computational Fluid Dynamics, elektromagnetism, elektriska maskiner,

PMSM, rotor, stator, lindningar, termisk, LPTN, körcykelanalys
Acknowledgements
I would firstly like to thank my supervisors at Chalmers, Alexey Vdovin and Emma
Grunditz for providing all the help and support one could ask for. A very special
gratitude to my supervisors at Volvo Car Corporation, Sriram Mandayam, Andreas
Andreasson and Raik Orbay for always encouraging me to push my limits. I have
learnt so much in this thesis, both academically and professionally thanks to all my
supervisors. A special mention to the manager of the team at VCC I was involved
with, Peter Berggren, for being a supporting overseer and for making me feel like a
part of his team. I would also like to thank Torbjörn Thiringer from Chalmers for
his guidance and feedback.

I am grateful to my parents and sister for putting up with my antics for all these
years, and also for giving me the opportunity to pursue my higher education. I will
always appreciate the emotional, mental and financial aid that they have provided
me. A special mention to my whole family and relatives for all their prayers and
wishes too.

I would like to thank my flatmates Adarsh, Akshay and Karthik, who have had
to deal with my incredibly bad puns for these two years and for many more to come!
Their friendship means the world to me. A big shout out to my friends and well
wishers too for always being there for me. A special bit of gratitude to my friend
Erik for helping me with the Swedish version of the Abstract.

Finally, I would like to thank the Person up above for all the opportunities that
have been bestowed upon me in this wonderful life.

Anthony Jayanath Vivek, Gothenburg, August 2019

vii
 Contents

Notations & Abbreviations xiii

List of Figures xv

List of Tables xix

1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Aim and Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Planning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4.2 Basic Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4.3 Design and Modelling . . . . . . . . . . . . . . . . . . . . . . . 5
1.4.4 Computational Fluid Dynamics (CFD) . . . . . . . . . . . . . 5
1.4.5 Lumped Parameter Thermal Network(LPTN or LPN) Model
Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4.6 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 Theory 7
2.1 Electric Motors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.1 Types of Electric Motors . . . . . . . . . . . . . . . . . . . . . 8
2.1.2 Modelling of PMSM . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Lumped Parameter Thermal Network . . . . . . . . . . . . . . . . . . 9
2.3 Heat Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3.1 Conduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3.2 Convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.3 Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.4 Computational Fluid Dynamics (CFD) . . . . . . . . . . . . . . . . . 14
2.4.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.2 Turbulence Modelling . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.3 Near Wall Modelling . . . . . . . . . . . . . . . . . . . . . . . 17
2.4.4 Rotation Motion . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4.5 Correlation of Results . . . . . . . . . . . . . . . . . . . . . . 20
2.4.6 Conjugate Heat Transfer (CHT) . . . . . . . . . . . . . . . . . 20

ix
Contents

3 Methods 23
3.1 CAD Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2 Rotational Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3 Mesh Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . 27
3.4 Basic EM Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.5 Thermal Resistances . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.6 Contact Resistances . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.7 Material Properties as a Function of Temperature . . . . . . . . . . . 33
3.8 Air-Gap Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.9 Solids Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.10 Active Winding Modelling . . . . . . . . . . . . . . . . . . . . . . . . 36
3.10.1 Windings Models . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.10.2 Anisotropic Material . . . . . . . . . . . . . . . . . . . . . . . 37
3.10.3 Windings CFD Modelling . . . . . . . . . . . . . . . . . . . . 38
3.11 End Winding Modelling . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.12 Cooling Jacket & Coolant Modelling . . . . . . . . . . . . . . . . . . 41
3.13 Complete Model & Validation . . . . . . . . . . . . . . . . . . . . . . 42
3.14 1D LPN Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.15 2D Thermal Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4 Results 51
4.1 Mesh Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . 51
4.2 Basic EM Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.3 Thermal Resistances . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4 Contact Resistances . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.5 Air-Gap Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.6 Active Windings Modelling . . . . . . . . . . . . . . . . . . . . . . . . 58
4.7 End Windings Modelling . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.8 Cooling Jacket & Coolant Modelling . . . . . . . . . . . . . . . . . . 61
4.9 Complete Model & Validation . . . . . . . . . . . . . . . . . . . . . . 63
4.10 1D LPN Model & 2D Thermal Model . . . . . . . . . . . . . . . . . . 65

5 Conclusions & Future Work 67

5.1 Component Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.2 Complete 3D Modelling . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.3 1D & 2D modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.4 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Bibliography 71

A Appendix A - Additional Images I

A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I
A.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II
A.2.1 Material Properties as a Function of Temperature . . . . . . . II

B Appendix B - Additional Information VII

B.1 Contact Resistances . . . . . . . . . . . . . . . . . . . . . . . . . . . . VII

x
 Contents

B.2 Air-Gap Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . VIII

B.3 Solids Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IX
B.4 Anisotropic Material . . . . . . . . . . . . . . . . . . . . . . . . . . . IX

C Appendix C - Physics Models XI

C.1 Air and Solids Physics . . . . . . . . . . . . . . . . . . . . . . . . . . XI
C.2 Coolant & Cooling Jacket Modelling . . . . . . . . . . . . . . . . . . XI

xi
Contents

xii
 Notations & Abbreviations

Notation Description Units

Tamb Ambient temperature K
h Convective heat transfer coefficient W/m2 − K
ρ Density kg/m3
µ Dynamic viscosity Pa − s
Tref Fluid reference temperature K
T0 Initial temperature K
ν Kinematic viscosity m2 /s
cp Specific heat capacity J/kg − K
σ Stefan-Boltzmann constant W/m2 − K 4
k Thermal Conductivity W/m − K
Rth Thermal resistance m2 − K/W

AC Alternating Current HTC Heat Transfer Coefficient

AW Active Windings HV High Voltage
CAD Computer Aided Design LES Large Eddy Simulation
CFD Computational Fluid Dynamics LH Left Hand
CHT Conjugate Heat Transfer LPN Lumped Parameter Network
DC Direct Current MRF Moving Reference Frame
DES Detached Eddy Simulation MW Moving Wall
DNS Direct Numerical Simulation PMSM Permanent Magnet Synchronous Motor
DOE Design of Experiments RANS Reynold’s Averaged Navier Stokes
EM Electric Machine RPM Rotations Per Minute
ERAD Electric Rear Axle Drive RSM Reynold’s Stress Model
EW End Windings SM Sliding Mesh
FEM Finite Element Method SST Shear Stress Transport
FVM Finite Volume Method TCR Thermal Contact Resistance
HMT Heat and Mass Transfer VCC Volvo Car Corporation

xiii
Contents

xiv
 List of Figures

1.1 ERAD motor[1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.1 Parts of a Permanent Magnet machine with an isometric view (left

and a cross-section view across the axial direction (right) [2]. . . . . . 8
2.2 The various thermal heat transfer paths in an electric machine[5]. . . 11
2.3 Heat transfer by conduction[20]. . . . . . . . . . . . . . . . . . . . . . 12
2.4 Thermal boundary layer in the fluid flowing over a heated surface[20]. 13
2.5 Velocity profile of a turbulent boundary layer [10]. . . . . . . . . . . . 17
2.6 An example of a prism layer with the core mesh in StarCCM+ [10]. . 18
2.7 The boundary interfaces between a fluid (left) and a solid (right) in
StarCCM+[10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.1 CAD model of the ERAD motor, created for this thesis. . . . . . . . 23
3.2 Velocity distribution in the rotor holes using a moving wall approach. 25
3.3 Velocity distribution in the rotor holes using a moving wall approach
(left) vs. a similar simulation with an MRF approach (right). . . . . . 26
3.4 The increase in mesh quality for the fluid domain from 0.1 million
(left) to 0.7 million (centre) to 1.9 million (right) cells. . . . . . . . . 27
3.5 The refinements made for the polyhedral mesher (left) and for the
trim mesher (right) for the simple conduction case using a rod. . . . . 28
3.6 The basic EM geometry with some dimensions. . . . . . . . . . . . . 29
3.7 The basic EM geometry with the fluid domain. . . . . . . . . . . . . . 29
3.8 The difference in mesh quantity without a glue region (left) and with
a glue region (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.9 The varying contact areas in a typical junction[25]. . . . . . . . . . . 31
3.10 Taylor vortices as seen in the axial cross-section of an EM (left)[28]
and uneven temperature distribution over the stator surface due to
Taylor vortices[30]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.11 The mesh in the air-gap with the rotor and stator also shown. . . . . 35
3.12 Windings model from the baseline model[2] and actual windings from
the ERAD machine. . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.13 The different windings models used in this thesis with the legend on
top. L-R: Complex; Baseline; Discretized; Anisotropic. . . . . . . . . 37
3.14 The thermal resistance configurations of different components in the
EM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.15 The boundary conditions for windings (left) along with a similar case
with a thermal resistance (right). . . . . . . . . . . . . . . . . . . . . 39

xv
List of Figures

3.16 The boundary conditions for EW (left) along with a similar case with
a thermal resistance (right). . . . . . . . . . . . . . . . . . . . . . . . 40
3.17 The cooling plate modelled for this thesis. . . . . . . . . . . . . . . . 41
3.18 Loss maps of copper losses (left) and rotor losses (right)[3][2]. . . . . 43
3.19 Different regions of the end windings and active windings in this thesis. 43
3.20 Full ERAD machine volume mesh in the XZ section plane (left) and
the YZ section plane (right). . . . . . . . . . . . . . . . . . . . . . . . 45
3.21 The 9-node LPN network from [3]. . . . . . . . . . . . . . . . . . . . 46
3.22 A scalar plot of the temperature distribution (left) and a vector plot
of the heat flux (right) obtained from the 2D rotor. . . . . . . . . . . 47
3.23 The 11-node LPN model created in GT-Suite for a 1/8th section of
the rotor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.24 The grid for a 1/16th section of the rotor with the different regions. . 49

4.1 The mesh sensitivity analysis for the fluids domain for the number of
cells in the domain vs. the temperature. . . . . . . . . . . . . . . . . 51
4.2 The temperature distribution over the rod for different mesh settings
and types. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.3 Specified y + HTC for the baseline model (left)[2] and the basic EM
model (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.4 Temperature distribution for the baseline model (left)[2] and the basic
EM model (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.5 Temperature distribution in the rotor with glue around the magnets
(left) and with a thermal resistance (right). . . . . . . . . . . . . . . . 54
4.6 Results of the DOE with varying R3 and R7 values and different
gearbox losses inputs. . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.7 Temperature distribution in the end windings for change in R3 and
R7 for a fixed gearbox loss. . . . . . . . . . . . . . . . . . . . . . . . . 55
4.8 Velocity vector in the air-gap annulus for 12000 RPM rotation rate,
when looking at an axial cross-section of the stator and the rotor. . . 56
4.9 Velocity vector in the air-gap annulus for 200 RPM rotation rate,
when looking at an axial cross-section of the stator and the rotor. . . 56
4.10 Temperature distribution over an arbitrary midway surface over the
air-gap annulus for 12000 RPM. . . . . . . . . . . . . . . . . . . . . . 57
4.11 Temperature distribution over a line probe over the air-gap annulus
for 12000 RPM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.12 Velocity vector in the air-gap annulus for 12000 RPM rotation rate,
when looking at an axial cross-section of the stator and the rotor. . . 58
4.13 Temperature distribution over a line probe over the air-gap annulus
for 12000 RPM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.14 Temperature distribution in the different models with similar bound-
ary conditions for active windings. L-R: Baseline model; anisotropic
model; discretized model. . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.15 Temperature distribution in the different models with similar bound-
ary conditions for active windings with a thermal resistance. L-R:
Baseline model; anisotropic model; discretized model. . . . . . . . . . 59

xvi
 List of Figures

4.16 Temperature distribution in the different models with similar bound-

ary conditions for active windings with a specific temperature im-
posed at one end of the windings. L-R: Baseline model; anisotropic
model; discretized model. . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.17 Temperature distribution in the different models with similar bound-
ary conditions for end windings. L-R: Baseline model; anisotropic
model; discretized model. . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.18 Temperature distribution in the different models with similar bound-
ary conditions for end windings with thermal resistances. L-R: Base-
line model; anisotropic model; discretized model. . . . . . . . . . . . . 60
4.19 Velocity vector in the fluid domain of the baseline model without
(left) and with thermal resistances (right). . . . . . . . . . . . . . . . 61
4.20 The comparison of the temperatures in different sensors for literature[36]
and the two CFD simulation cases. . . . . . . . . . . . . . . . . . . . 62
4.21 The temperature distribution in cooling plate. . . . . . . . . . . . . . 62
4.22 The temperature distribution in full ERAD motor for the rig operat-
ing conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.23 The temperature distribution in full ERAD motor for the second case. 64
4.24 The temperature distribution in the 1/16th section of the rotor from
the 2D model with the colour bar showing a representative tempera-
ture scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

A.1 The initial planning Gantt chart. . . . . . . . . . . . . . . . . . . . . I

A.2 The thermal conductivity of copper varying with temperature. . . . . II
A.3 The thermal conductivity of lamination steel as a function of temper-
ature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II
A.4 The thermal conductivity of the shaft steel varying with temperature. III
A.5 The variation of dynamic viscosity with temperature of the coolant. . IV
A.6 The variation of thermal conductivity with temperature of the coolant. IV
A.7 The variation of density with temperature of the coolant. . . . . . . . V
A.8 The variation of dynamic viscosity with temperature of the air. . . . . V
A.9 The variation of density with temperature of the air. . . . . . . . . . V

xvii
List of Figures

xviii
 List of Tables

3.1 Differences in dimensions from actual ERAD motor and Anil’s geometry. 24
3.2 Operating parameters for sensitivity analysis. . . . . . . . . . . . . . 27
3.3 Operating parameters and heat loss input conditions for the simple
case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.4 Operating Parameters for sensitivity analysis. . . . . . . . . . . . . . 33
3.5 Boundary conditions for CFD simulation of literature cooling plate. . 42
3.6 The volume fractions of different winding regions. . . . . . . . . . . . 44
3.7 Operating parameters for the rig setup comparison simulation. . . . . 44
3.8 Volume cell count in different regions of the simulation. . . . . . . . . 45
3.9 Operating parameters from the baseline simulation [2]. . . . . . . . . 46
3.10 Calculated Biot number for the ERAD motor. . . . . . . . . . . . . . 47

4.1 Differences in results for temperature and Specified y+ HTC of the

basic EM model from the baseline model. . . . . . . . . . . . . . . . . 53
4.2 Differences in results for temperature across the two simulations. . . . 53
4.3 Results from CFD compared against the rig measurements for similar
operating conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.4 Temperatures in different regions of the ERAD motor compared be-
tween the baseline model[2] and the simulation from this thesis. . . . 64
4.5 HTC in different regions of the ERAD motor compared between the
baseline model and the simulation from this thesis along with the
reference temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.6 Temperature in different regions from a steady state simulation of the
1D LPN model of the 1/8th rotor. . . . . . . . . . . . . . . . . . . . . 65

B.1 TCR calculation assumptions and results for rotor and shaft . . . . . VII
B.2 TCR calculation assumptions and results for stator and cooling jacket VIII
B.3 Material properties for solids in the ERAD machine . . . . . . . . . . IX

C.1 Physics Continuum for the Basic EM Case . . . . . . . . . . . . . . . XI

C.2 Physics Continuum for the Coolant & Cooling Jacket Modelling Case XI

xix
List of Tables

xx
 1
Introduction

1.1 Background
With the global automotive market moving towards electric mobility, highly effi-
cient electric drive systems are required for sustainability of the environment and
for individual automotive manufacturers. With a higher efficiency of the drivetrain
system, a greater mileage can be extracted for the vehicle. For a given unit of energy,
with greater range coverage of the vehicle, a reduced carbon footprint of the vehicle
is achieved. Such vehicles can also be regarded with respect in industry, thereby
improving the manufacturer’s market credibility.

Some of the major components in an electric drive system include accumulators,

inverters, motors, transmission, etc. A way to determine the efficiency of the com-
plete system is by running drive cycle analysis on the system. For this, different
losses in the system need to be determined. Of the losses in the electric system, the
mechanical losses are one of the larger losses. Most of these mechanical losses are
transformed into thermal losses. Armature losses due to current in the windings,
and ohmic losses due to the resistance in the windings, also lead to thermal losses.
Therefore, by reducing the thermal losses in the system, the efficiency of the overall
system can be improved. To do this, good thermal path predictions are required to
design and validate cooling methods for those parts. This thesis solely concentrates
on the electric motor, its thermal paths and methods involved in cooling the motor.

Thermal paths can be simulated with complex geometries and reasonably real op-
erating conditions using 3 dimensional (3D) Computation Fluid Dynamics (CFD)
simulations. These are computationally very expensive, but with proficient mod-
elling of the physics, they provide quite accurate results. Performing complex drive
cycle analysis with transient conditions and variable operating parameters will take
a large amount of time and resource. Hence, the machine can be represented as a
one dimensional (1D) or two dimensional (2D) model taken from literature, mea-
sured data or CFD results.

1D models are very fast and do not require high computing power, but on the
other hand are as reliable as the formulations used. The motor can be modelled
using a Lumped Parameter Network (LPN) approach, for the 1D modelling. LPN
models consider conducting solids to be lumps with thermal resistances between
different lumps. In addition to conduction, convection (both these phenomenon will

1
1. Introduction

be explained in the Theory chapter) can also modelled as thermal resistances as a

function of the properties of the fluid or the operating parameters. Thermal ca-
pacitances are modelled for performing transient simulations. The motor model is
simulated with 3D CFD software StarCCM+, and these results are used to improve
the existing 1D model.

Accurate and efficient (in terms of performance and resources used) modelling of
the thermal paths in an electric motor is a vital component in electric motor inter-
nal thermal analysis. In this thesis, the thermal model of the motor will be modelled
in steps of increasing complexity and also as a modular component. Firstly, basic
simulations will be carried out with a concentric parts motor concept. This way, the
simulations act as a tutorial in this thesis, for thermal simulations in StarCCM+.
This is followed by modelling separate regions of the motor in conjunction with
empirical data or experimental data from literature. Finally, all these models are
integrated together to run complete motor simulations. The motor in question is the
Electric Rear Axle Drive (ERAD) motor used in the Volvo XC90 cars. The CAD
is modelled based on dimensions obtained from a disassembled motor, at Chalmers.
Before it was dismantled, the motor was run in a rig at Chalmers, to obtain temper-
ature readings on the end windings, at certain operating conditions. The complete
CFD model created for the motor will be compared against the rig measurements,
by simulating at similar operating conditions.

Figure 1.1: ERAD motor[1].

Good corelations of the CFD data with the rig measurements will provide an accu-
rate base to calibrate or validate the LPN model. A well modelled LPN network will
provide correct thermal resistance paths and hence can be utilized for drive cycle
analysis with reliable results. This better modelling of the LPN model not only helps
industrial technical data analysis for various motors used in the automotive sector,
but also from a technical research point of view, to understand the advantages and
shortcomings of both the 3D CFD simulations and the 1D LPN Networks. This

2
 1. Introduction

understanding can help recognize where the LPN can be effectively used, and where
it falls short in terms of accuracy of results.

1.2 Previous Work

As stated before, this thesis is a continuation thesis from [2]. It primarily involved
3D CFD modelling for the heat transfer in electric machines. The thesis work was
carried out under the supervision of Chalmers and Volvo Car Corporation, same as
in case of this thesis. The motor in question is the same Volvo ERAD motor. In
[2], the modelling was done first with a simple model, and then with a 1/8th portion
with symmetry conditions imposed and finally on the complex model. Specific areas
such as end-windings or air-gap were not considered for detailed modelling and a
general modelling was carried out. The simulations were run for various operating
parameters and the obtained results were provided to the LPN model which was
developed from [3] in a project at VCC. The coolant flow had been considered to be
laminar for all flow rates in [2].

The LPN model which is being improved in this thesis is based on [3]. These LPN
models implement empirical formulations for different regions of the motor based on
literature. Other parameters or input data are obtained from 3D CFD modelling
conducted by various students or researchers. The LPN model which this thesis is
aiming to improve has a simple formulation for the end-windings region and does not
have axial nodes (completely radial 1D model). Hence, this thesis aims to improve
these regions in 3D CFD and then implement those results appropriately into the
1D LPN model.

1.3 Aim and Scope

The aim of this thesis is to improve the modelling of the LPN of the electrical ma-
chine. By using the results of the 3D CFD simulations, the LPN can be calibrated
or modified to provide results more relatable to existing measured data. An alter-
native method is to model the LPN using new approaches. This way, the LPN will
be independent of test data or track measurements as possible. This new method
will implement 2D thermal models or simple empirical calculations to complement
the LPN model and use minimal 3D CFD computations.

The scope of the thesis includes:

• Pre-study - Literature review, journals and previous theses

• Understanding the software - Basic CHT and CFD simulations, mesh sensitiv-
ity analysis, physics settings and rotational simulations
• Basic EM model - Basic representative EM model consisting of concentric
cylinders
• Air-gap study - Modelling of air-gap in CFD and comparison to literature

3
1. Introduction

• Windings modelling - Modelling the active and end windings using multiple
approaches and selecting suitable model for the thesis
• Coolant modelling - Modelling of cooling jacket and coolant and comparing
to literature; Modelling materials with physical properties as a function of
temperature
• Thermal contact resistance - Estimating contact resistances in EM and im-
plementing in CFD; Replacing complex or small geometry with equivalent
thermal resistances
• Full EM simulations - ERAD machine 3D CFD simulations and comparison
of results to rig measurements and previous theses
• Basic LPN modelling - Simple LPN models using GT-Suite based on 2D sim-
ulations

1.4 Planning
The thesis is split into various phases and each phase is elaborated in their respective
sections in the Methods Chapter.

• Pre-Study and Literature Review

• Electromagnetic Losses Modelling
• Basic Model
• Design and Modelling
• EM Component Modelling and Validation
• Full EM Simulations and Validation
• LPN Modelling and Calibration

The thesis is planned for a period of 5 months, stretching from Week 3, 2019 to
Week 25. A detailed planning chart or a Gantt chart is prepared describing the
methods or processes in the thesis, and their expected timeline. This Gantt chart
can be seen in Figure A.1 in Appendix A. Halfway through in the thesis, the plan-
ning document has been updated to make suitable changes. The plan to possibly
couple 3D CFD and 1D LPN models is eliminated. Instead, a bottom-up approach
has been incorporated in the thesis wherein each component of the motor has been
modelled individually. These models are validated against thermodynamic calcula-
tions, empirical data or literature. The same physics, mesh setting and modelling
approach can be utilized for the components in the complete machine modelling.
Along with these components, other phenomenon have also been modelled at the
same time, before incorporating it in the full machine simulations. These include
thermal resistance calculations to replace small, complex geometries, rotational sim-
ulations, mesh sensitivity analysis and adding material properties as a temperature
dependant function. All of these will be defined in the Methods chapter.

1.4.1 Literature Review

This phase involves going through the theory involved in this thesis. This includes
textbooks, articles, journals, previous thesis reports (Master thesis and PhD. re-

4
 1. Introduction

ports) and online content. Apart from the literature review, a pre-study phase is
conducted involving planning the thesis, defining a scope and targets, listing soft-
ware to be used in the project, and resource gathering (software licenses and access
at VCC and Chalmers, documents, supercomputer access, desk allocation at VCC
and Chalmers).

1.4.2 Basic Modelling

This phase involves going through software tutorials, running mesh sensitivity anal-
ysis, choosing physics approaches for future simulations and finally using the knowl-
edge from the above listed items in running a simulation for a basic EM model
involving concentric cylinders representing the rotor yoke, rotor magnets, active
windings and stator yoke.

1.4.3 Design and Modelling

Here, the geometric modelling of the machine and other components for component
level simulations are looked into. Most of the design is carried out based on the
ERAD machine, available at Chalmers. The dimensions are hand measured and
then reproduced in a CAD software and further cleaned in ANSA.

1.4.4 Computational Fluid Dynamics (CFD)

In this stage, the computational simulations for different components or for multiple
components are conducted. These simulations also include Conjugate Heat Transfer
(CHT) to involves heat transfer between solids and fluids. They are mainly con-
ducted in StarCCM+ and sometimes in other CFD software such as OpenFOAM.

1.4.5 Lumped Parameter Thermal Network(LPTN or LPN)

Model Analysis
Once CFD simulations are conducted, the obtained results can be used to assess
existing LPN models or to create new models for the electric machine. These models
are processed using MATLAB or GT-Suite.

1.4.6 Validation
Every result obtained in this thesis needs to be validated against a reliable source.
This source may be a theoretical formulation, an empirical solution provided by con-
ducting experiments, or from real life data such as rig or drive-cycle measurements.
Validation approves the physics models used in simulations and can further be used
to model similar components.

5
1. Introduction

6
 2
Theory

In this section, the various components, modelling techniques and the physics behind
these methodologies are elaborated. Electric machines, different types of machines
and their various components are looked into here. The physics behind heat transfer
and the modelling techniques used in this thesis (LPN and CFD) and the theory
behind the modelling techniques are also explained in this chapter.

2.1 Electric Motors

An electric motor is a machine that can convert electrical energy to mechanical
energy or vice versa. This is done by the interaction of magnetic forces and cur-
rent passing through conductors to produce mechanical forces. The motor consists
usually of a rotor and a stator, with either electrically conductive windings in both
components or windings in one and a permanent magnet in the other. The stator
is fixed in nature, usually to an external casing, and the rotor is placed on a shaft
with bearings, that transmits the mechanical forces as a rotational motion[9]. The
forces are created based on Fleming’s left hand (LH) rule which states that if the
current flows in the direction of the middle finger, the magnetic field is created in
the direction of the index finger, and the motion in the direction of the thumb.

Figure 2.1 shows the motor and the various components which are listed below.
• Shaft - It transfers mechanical force from rotor to other components such as
transmission or pulleys.
• Rotor - It houses the conductive windings or magnets on which the (elec-
tro)magnetic forces act to produce mechanical forces. It is made of electro-
magnetic laminations separated by an electrically insulating bonding material.
• Magnets/Windings - They are the elements on which the electromagnetic
forces from the active windings act to create mechanical forces. These ele-
ments can be permanent magnets housed in the rotor or they can be in the
form of excited windings which produce a similar magnetic field.
• Air Gap - Between the rotor cylinder and the stator cylinder, there exists a
small gap, occupied by air.
• Stator - Similar to the rotor, it houses the conductive active windings which
produce the electromagnetic field. It is made of electromagnetic laminations
separated by an electrically insulating bonding material.
• Slot - The slots in the stator are the housing channels for the active windings.

7
2. Theory

• Slot Liner - The slot liners are electrically insulating lining material placed
between the stator and the windings.
• Active Windings - They are windings of conductive material (usually copper)
through which current is passed. They are a part of the windings present in
the slots in the stator and are wound around to neighbouring slots.
• End Windings - They are the part of the windings that come out of the sta-
tor and go around back into neighbouring slots. They do not participate in
producing any torque in the machine and hence are not considered as active
windings.
• Cooling Channel - Due to high thermal loads in the electric machine, a cooling
jacket or channel is run around the stator and sometimes around the windings
too. They carry a coolant fluid around to extract heat from the stator yoke.

Figure 2.1: Parts of a Permanent Magnet machine with an isometric view (left
and a cross-section view across the axial direction (right) [2].

2.1.1 Types of Electric Motors

The most common electric motor in the automotive industry is the AC motors. In
these motors, the windings are activated by alternating current supply (single phase
or multi phase). By using alternating current, the need for switching mechanisms
is eliminated. Some examples of AC motors are Induction Motors and Synchronous
Motors - These are motors where the shaft rotation is synchronized with the supplied
current frequency. Some of the common synchronous motors include reluctance
motors, inductance motors and the one used in this thesis - Permanent Magnet
Synchronous Motor (PMSM).

2.1.2 Modelling of PMSM

The PMSM motor (or any motor) design is conceived based on the following aspects:

• Performance Characteristics
• Duty Cycles

8
 2. Theory

• Service Factors
• Motor Efficiency
• Thermal Aspects
• Packaging
• NVH (Noise Vibration and Harshness)

Using these parameters, a motor can be designed in multiphysics software which

can be coupled with electromagnetic, thermal and mechanical solvers to simultane-
ously design an optimum motor for all the above stated parameters [6]. The motor
can also be represented as Lumped Parameter Thermal Networks (LPTN or LPN),
which will be elaborated in the next section.

For a motor, the electrical power, the loss PC u , which represents the windings losses,
is estimated by the equation 2.1.

PC u = I 2 ∗ RC u (2.1)

Where,
I is the current flowing through the cable, A
RC u is the resistance in the copper cable given by RC u = ρ ∗ l/A
ρ is the density of the cable in kg/m3
l and A are the length and the cross-sectional area of the cable

2.2 Lumped Parameter Thermal Network

LPN is a thermal modelling approach where solid objects are lumped into single
blocks for simplified thermal calculations, wherein the Biot number (Bi) provides
the validity of the LPN. Biot number is a dimensionless number which represents
the ratio of conductive heat resistance within a solid to convective heat resistance
at the solid’s boundary. Biot number is given by the equation 2.2.
h ∗ Lc
Bi = (2.2)
k

Where,
h is the heat transfer coefficient of the fluid at the surface of the solid, W/m2 − K
Lc is the characteristic length of the solid, m
k is the thermal conductivity of the solid, W/m − K

For an LPN to have reasonable results, the Bi < 0.1, as in this case [7][8], the
temperature completely within the solid will be quite uniform, as compared to the
temperature at the surface of the solid. With a higher Bi value for different compo-
nents, it implies that the parts are being effectively cooled, but that in turn signifies
that the temperature distribution in the part is not similar throughout and hence,
the LPN for that component needs discretization[8].

In the LPN used in this thesis, radiation is ignored as the temperature difference

9
2. Theory

(raised to the power of 4) between adjacent components are negligible as compared

to the Stefan-Boltzmann constant (5.67 ∗ 10−8 W/m2 K 4 ) used for radiation emission
calculations as per Equation 2.7 [3][4]. Another reason to ignore the effects of radi-
ation are because the  of the materials in the ERAD have very low values.

These lumps are connected by nodes, placed usually at the centre of each lump.
These connections created are then modelled as thermal resistances and capaci-
tances based on the thermal properties of the material of the connected lump. All
capacitances are based on the material specific heat values and are used for transient
calculations. For solid-solid connections, the thermal resistances are based on the
conductive resistance between the two solids in question.

For solid-fluid connections, the thermal resistances are based on the convective ther-
mal resistance of the fluids. These convective resistances are a function of heat
transfer coefficients between the solid and the fluid, which in turn are a function of
various other parameters (material properties, fluid velocity, etc.). In the case of
the electric machine, for the air-gap for instance, the heat transfer coefficient is a
function of air-gap thickness, angular velocity of the rotor, the air properties and the
rotor dimensions[28]. Similarly, for the coolant jacket, the heat transfer coefficient
in this case is a function of fluid flow rate, fluid properties, operating speed and loads.

LPN for a motor is created, on literature for different components using empirical
formulations or thermal formulations for conduction and convection. Some of the
parameters cannot be determined by thermal calculations or empirical formulations
and hence require CFD simulations. The heat transfer around the end windings for
example, require CFD simulations for various load conditions and operating speeds.
These values are tabulated and extrapolated for the whole range of loads and oper-
ating speeds.

Once this LPN for a motor is modelled, in simulation software such as MAT-
LAB/Simulink or GT-Suite, it can be added to system models which include bat-
teries, inverters, transmission and vehicle dynamics. Using this model, drive cycle
simulations can be run. Drive-cycle analysis is used to virtually simulate the condi-
tions of a car based on different conditions such as highway driving, urban driving,
high performance driving, range calculations, etc. Hence, a well calibrated LPN
can be matched quite closely to real life drive cycle data, and this is very useful in
testing new components, or even in designing various components of an electrical
propulsion system.

2.3 Heat Transfer

Heat transfer between two bodies is the rate of change of energy between the bodies
due to an existing temperature gradient between them. Here, the laws of thermo-
dynamics dictate the heat transfer. The first law states that for two closed systems,
the rate of energy entering a system must be equal to the rate of energy disposed

10
 2. Theory

by the other system. The second law states that heat can only be transferred due
to a temperature gradient. Hence, larger the temperature gradient, greater is the
heat transfer[20].

For a unit mass of any given substance, the energy required to increase it’s tempera-
ture by 1 degree is called it’s specific heat. There is specific heat at constant pressure
and specific heat at constant volume represented by cp and cv respectively. Usually,
for low Mach numbers, fluids are usually considered to be incompressible. Similarly,
most solids at working temperatures are considered to be incompressible. For solids
and fluids, cv and cp are the same[4]. Then the internal energy of materials, U can
be defined as:

∆U = mc∆T (2.3)

Where,
m is the mass of the system in kg,
c is the specific heat of the material c = cv = cp J/kg-K
∆T is the change in temperature in K

Figure 2.2: The various thermal heat transfer paths in an electric machine[5].

Heat energy can be produced by converting different forms of energy. Electrical re-
sistance in the coils can be converted to heat. The friction between air and rotating
components, however minute, can create heat. Varying electromagnetic fields in-
duced in components can create heat. Hence, it is vital to estimate the heat loss for

11
2. Theory

different components in an electric machine. Improper extraction of the generated

heat in the machine may lead to unwanted consequences, such as loss of performance
and/or damage of components. To avoid this, proper thermal paths must be ap-
proximated in the machine and hence, appropriate cooling solutions must be devised
based on the paths. The Figure 2.2 shows the thermal transfer paths in an electric
machine[5]. The red arrows indicate conduction between solids, the green arrows
indicate convection in the fluids and the purple arrows represent radiation in the
components. These are the three heat transfer mechanisms and will be explained in
the subsections below.

2.3.1 Conduction

Figure 2.3: Heat transfer by conduction[20].

Conduction is a heat transfer phenomenon which occurs at the molecular level and
occurs due to energy transfer from molecules having greater thermal energy to ones
having lesser thermal energy[20]. The ability of a material to transfer energy through
conduction is represented by a coefficient called thermal conductivity or k having the
units W/m − K. Consider a solid with a temperature gradient between the walls, as
shown in Figure 2.3. The 1D heat flux, qx ” in W/m2 , through a material of thickness
L m, due to conduction, for a steady-state heat transfer is given by the equation 2.4.

−k
qx ” = (T2 − T1 ) (2.4)
L

The conductivity is negative because the heat is transferred in the direction of

lower temperature. For a time-dependant heat transfer calculation in solids, the
heat transfer equation is defined with respect to the thermal diffusivity, α, of the
material. The heat equation is given by:

12
 2. Theory

dT
= α∇2 T (2.5)
dt

Where,
Thermal diffusivity, α = ρckp
ρ is the density of the material kg/m3
cp is the specific heat at constant pressure for the material, J/kg − K

2.3.2 Convection

Figure 2.4: Thermal boundary layer in the fluid flowing over a heated surface[20].

For fluids in contact with solids, the heat is transferred by two methods. There is
conduction at the walls between the solid and the fluid interface[4]. The majority of
the heat is then transferred due to the bulk motion of the fluid. Consider a flowing
fluid on a heated surface with a free stream velocity temperature of of u∞ and T∞
respectively, as seen in Figure 2.4.

At the walls, the fluid velocity will be zero, due to a no-slip condition. Going
further into the free stream, the velocity normalizes to u∞ . This phenomenon is
visualized as velocity boundary layer. Similarly, due to the surface temperature Ts
at the walls, the fluid temperature at the walls will be close to Ts when the flow
conditions are time-independent (steady state). Moving away from the walls, the
temperature profile will be similar to an inverted velocity boundary layer and will
normalize to T∞ due to bulk movement of fluid, or convection. This boundary layer
is known as the thermal boundary layer[20].

Convection can be of two types: forced or natural. When the fluid motion is due
to an external source, such as a fan or a pump, the convection is known as forced

13
2. Theory

convection. When the temperature changes in a fluid, according to the gas law
P V = mRT , the density varies as well. This variation in density causes buoyant
forces to act on the fluid, inducing motion. This induced motion causes a convection
current in the fluid and this type of convection is called natural convection.

The convection heat transfer is defined by the equation 2.6.

q” = h(Ts − T∞ ) (2.6)

Where,
h is the heat transfer coefficient in W/m2 − K and is defined by the type of flow
(forced or natural) and by the properties of the fluid.
Ts is the surface temperature at the boundary in K
T∞ is the ambient fluid temperature in K

2.3.3 Radiation
Any substance having a temperature greater than 0 K emits thermal radiation in
the form of electromagnetic waves. Radiation does not need a medium for heat
transfer, contrary to conduction or convection, in which the medium of propagation
of heat is vital[20]. The heat flux E, emitted due to radiation from a body having a
surface temperature Ts is:

E = σTs4 (2.7)

Where,
 is the emissivity coefficient of the body ( = 1 for a black body)
σ is the Stefan-Boltzmann constant (σ = 5.67.10−8 W/m2 − K 4 )

In the electric machine, although a few hot surface temperatures (rotor, magnets,
end-windings) are present, the emissivity coefficients are quite low for the materials
(0.01-0.1) due to their manufacturing method[21]. Hence, the effects of radiation
can be ignored for simulations inside an electric machine. On the other hand, heat
transfer from the electric machine to the outside will need to take radiation into
consideration.

2.4 Computational Fluid Dynamics (CFD)

The fluid flow is governed by the Navier-Stokes equation. This equation is based
on the physical properties of the fluid. With advances in computational sciences,
the governing equations are modified and can then be solved numerically to esti-
mate the fluid flow. Computational Fluid Dynamics is the field of fluid dynamics
to approximate fluid flow using numerical codes. These codes may be commercial

14
 2. Theory

or in-house in nature. The commercial code provided by Siemens called StarCCM+

v1306[10] is utilized in this thesis.

The simulations are Eulerian in nature meaning that the fluid is observed in a
period of time in a given space. The computational software can solve for the en-
ergy equation as well and hence thermal calculations are conducted as well. For the
fluid modelling, various turbulence models are available and for for discretizing the
viscous layer at the boundary, many near wall modelling methods are present. To
simulate the effect of rotation or motion, various methods are present, explained in
detail below.

In an electric machine, CFD can be used extensively to analyze the complex fluid
flow of the air around the machine, and also the coolant flow. This flow cannot be
predicted by simple empirical models, and hence CFD is preferred. The CFD model
can then be used to create, validate or improve the LPN model of the machine.

2.4.1 Governing Equations

Fluid motion and heat transfer in materials are governed by the three laws of conser-
vation. Using these conservation equations, the Reynold’s-averaged Navier-Stokes
(RANS) equation has been formulated to define fluid motion. This equation is
named after the physicists, Claude-Louise Navier and George Stokes and is as seen
in Equation 2.9 for incompressible flows.

• Conservation of Mass - For incompressible fluids, mass flow into the system
is always equal to the mass flow out of the system as seen in the equation 2.8.

dρ
+ ρO.ui = 0 (2.8)
dt
Where,
ρ is the density of the fluid in kg/m3
t is the time in s
ui is the velocity component in the x,y and z directions in m/s

• Conservation of Momentum - For incompressible fluids, within an isolated

system, total momentum remains a constant. This is characterized in the
equation 2.9.
dui ∂p ∂τj i
ρ =− + + ρFi (2.9)
dt ∂xi ∂uj
Where,
Fi is an external forces in N
τj i is the shear stresses in the respective directions in N/m2
p is the pressure in Pa

• Conservation of Energy - In a system, the total amount of energy remains

constant and can only change form within the system. This is defined in the

15
2. Theory

equation 2.10.
de ∂ ∂T ∂ui ∂uj
ρ = ρq̇ + (k ) − p( + τij (2.10)
dt ∂xi ∂xi ∂xi xi
Where,
e is the internal energy of the system
q̇ is the heat flux entering or leaving the system
T is temperature in the system

By solving these equations, with some simplifications or modifications (for turbu-

lence characteristics), the properties of the flow can be obtained. The velocity and
temperature distributions in a fluid field can also be determined and thereby com-
plex aerodynamic and thermal problems can be solved.

2.4.2 Turbulence Modelling

Majority of fluid flow problems are turbulent in nature[12]. In an electric machine,
the air experiences a lot of turbulence due to the spinning rotor. The air-gap contains
Taylor vortices, which will be explained in Section 3.8, which causes entrainment at
the ends, thereby creating turbulence around the end windings as well. The coolant
flow is typically turbulent due to the sharp bends[13][14][15] and high mass flow
rates. To capture all these turbulent flows, different turbulence modelling physics
are required, based on the mesh conditions, wall conditions, accuracy and resource
availability. A property, for example velocity vi , in turbulent flows is usually split
into two parts for calculations: a time averaged part, vi and a fluctuating term, vi0 .

vi = vi + vi0 (2.11)

Some of the characteristics of a turbulent flow[16] are irregularity, diffusivity, 3 di-

mensional flow and dissipation. There are various turbulence models available such
as Reynold’s-averaged Navier Stokes (RANS) based models, Large eddy simulation
(LES), Detached eddy simulation (DES) and Direct numerical simulation (DNS) to
name a few. In this thesis, since the scope is limited to steady-state simulations,
the focus is primarily on the k- and k-ω RANS based models. They are both two
equation models and consider that the turbulence field is related to the mean veloc-
ity field using a linear eddy viscosity coefficient. This assumption is based on the
Boussinesq approximation[12]. The two equation models account for the turbulent
kinetic energy k and depending on the model, turbulent dissipation  or specific tur-
bulent dissipation rate ω. The specific models used in this thesis are the Realizable
k- model and the Shear Stress Transport (SST) k-ω model.

The Realizable k- model, a modification of the k- model, is a good model for free-
shear layer flow, where the flow does not have large adverse pressure gradients[17][18].
Hence it is suitable for large, open air simulations, but not so much for compressors

16
 2. Theory

and closed off geometries[19]. As recommended by the StarCCM+ manual, for the
Realizable k- model, the high-y + and all-y + wall treatment should provide reason-
able results.

The shear stress transport (SST) k-ω model combines the best parts of the k-ω
model and the SST model. The viscous sub-layer boundary layer formulation from
the k-ω model provides good results all the way down to the wall without the need
for further damping functions to eliminate the blocking effects of the wall. The SST
model can switch to a k- model in free stream approximations. The StarCCM+
recommended wall treatment for this model is the high-y + and all-y + wall treatment.

2.4.3 Near Wall Modelling

With any flow simulations around a body, proper prediction of flow around the
walls is paramount to the development of the flow almost everywhere else around
the body. At the wall, the velocity of the fluid relative to the wall is zero due to the
no-slip condition of the wall. Going further into the fluid stream, the velocity grad-
ually increases and then saturates close to the velocity of the fluid free stream[16].
The velocity profile of the fluid from the boundary layer can be seen in the image
below.

Figure 2.5: Velocity profile of a turbulent boundary layer [10].

• Log layer - This layer is dominated both by viscous and turbulent effects.
• Buffer layer - It is the transitional layer in between the log layer and the sub
layer.
• Viscous sub layer - In this layer, the fluid is dominated by viscous effects. The
flow is majorly laminar in nature.

The distance of the inner layer is measured by a non-dimensional wall distance called
y + . It is given by the equation 2.12.
u∗ ∗ y
y+ = (2.12)
ν
17
2. Theory

Where, q
u∗ is friction velocity at the nearest wall given by τρw in m/s
τw is the wall shear stress in N/m2
ρ is the density of the fluid in kg/m3
y is the nearest distance to the wall in m
ν is the kinematic viscosity of the fluid in m/s2
In the viscous sub layer, the law of the wall takes precedence and the velocity profile
can be estimated by the equation 2.13.
1
u+ = ln(y + ) + B (2.13)
k
Where,
u+ is dimensionless velocity
k is the von-Karman’s constant ≈ 0.41
B is a constant ≈ 0.51

Many CFD turbulence models usually do not require auxiliary treatment models in
the free stream volume, but at the viscous sub layer, where the law of the wall takes
precedence, additional wall treatment functions are necessary. For a fluid flow sim-
ulation, at the boundary layer, prism layers can be added to resolve the boundary
layer better, as the flow very close to the boundary layers are usually tangential or
laminar, as stated before. The number of prism layers and the total height of the
prism layers is usually specified based on the turbulence model and their recom-
mended wall treatment.

Figure 2.6: An example of a prism layer with the core mesh in StarCCM+ [10].

StarCCM+ has standard wall functions for the viscous sub layer and blended wall
functions which blends the log layer and the viscous sub layer to represent the
buffer layer[10]. In StarCCM+, 4 wall treatment methods are present, namely Low-
y + , High-y + , All-y + and 2 layer-y + wall treatment. The All-y + wall treatment is a
composite treatment approach that uses a High-y + wall treatment for coarse meshes

18
 2. Theory

(y + ≥ 30) and the Low-y + wall treatment for fine meshes (y + ≤ 5). Even if the
mesh falls within the buffer layer (5 ≤ y + ≥ 30), it is still formulated to give rea-
sonably good results. In this thesis, the All-y + wall treatment is mainly used for
the turbulence models as the mesh varies in quality depending on the component of
the electric machine (fine mesh in the air-gap and coarse mesh in the air around the
rotor shaft). Since the treatment approach gives good results for all y + values in the
boundary layer, it is a good complement to the turbulence models used in the thesis.

In CFD, wall treatments are available, to account for the boundary layer. To prop-
erly resolve the boundary layer, prism layer cells are used, in StarCCM+. This
mesher is used with a base volume mesher and adds prismatic layers on the wall
surfaces. The number of layers, growth rate and layer height can be defined in
StarCCM+ [10]. The documentation recommends the wall y + value required for
different turbulence models and their respective wall treatments selected. Based on
these recommendations, the number of prism layers can be selected and the prism
layer height is approximated based on Equation 2.12. Here, the value of y obtained
from the recommended y + value will be the first prism layer cell height from the
surface.

2.4.4 Rotation Motion

To simulate the rotational motion of the shaft and the rotor in the machine, there
are three ways in StarCCM+. They are:

1. Rotating wall
In this approach, a tangential velocity component is applied to the specified
walls. This approach is best suited for walls which are fully or nearly tan-
gential to the flow direction. In case of flat surfaces, a slip condition can be
applied to the walls which will not decelerate the flow near the boundary.

2. Moving Reference Frame (MRF)

Motion can be specified to a specific region in the domain, by applying a rel-
ative motion to the laboratory coordinate system. This motion can be based
on a relative reference coordinate system. The motion can be translating or
rotating or both in nature.

3. Moving Mesh
In this method, the mesh cells are translated or rotated with every time step
and updated to imitate the motion of the region on which it is applied. It
is very time and resource consuming as for every step, the mesh is updated
and numerous inner iterations are conducted for each time step. This sort of
approach is only for a transient simulation and hence is not considered in this
thesis.

19
2. Theory

2.4.5 Correlation of Results

With the use of various modelling techniques in a simulation environment, there will
be a lot of differences with the results for a given region of the electric machine. For
example, for the windings, multiple approaches are used to model them. Different
models may provide different results and to verify these results, different methods
of validation will be required. A good way is to analyze the equations behind the
modelled physics in the simulation and then compare them to empirical equations or
to theory. This outlook to the problem is called correlation of simulation to theory.
In this thesis, thermal resistances and simple conduction cases are correlated with
empirical formulation from theory.

Researchers over the years have tried and tested most approaches or modelling
techniques for electrical machines, and will have documented it in journals, confer-
ence papers or other literature. This database is a very good source for correlation
with simulations. Although the geometry or case may not completely match for
this thesis, it will usually be very similar and have the same dominating physics
conditions. Hence, modelling a test case similar to literature can be very useful to
assimilate mesh settings, boundary conditions, turbulence models and other physics
settings. A good correlation will imply that a similar modelling of the components
in this thesis will provide very good results.

Correlation of literature data to simulation results is used in this thesis for the
air-gap modelling, contact resistance approximation, LPN modelling and coolant
system modelling. All these cases will be explained in detail in the Methods chapter.

The most reliable validation of a modelling technique is by comparing the simu-

lation results directly with a real life setup of the simulation. This is usually con-
ducted in a rig setup or by emulating the conditions in reality. The ERAD motor
has been tested at Chalmers in a rig setup, and also drive cycle analysis from test
runs at the Hällered test track are available. This data can be compared directly
to simulations having similar boundary conditions. Usually this is very expensive
or time-consuming to setup and execute and is not suitable for concept validation
unless the concept is certain to provide good results.

2.4.6 Conjugate Heat Transfer (CHT)

Conjugate heat transfer is the solving of a thermal problem involving both solids
and fluids, by solving for the heat transfer at the interface of the two different ma-
terials. The temperature field in the solids is calculated by empirical conduction
equations and the temperature at the walls of the solids will then act as the bound-
ary temperature conditions for the fluids as well. The temperature field in the fluids
is calculated by the energy equations, as stated in Equation 2.10. In StarCCM+,
the heat transfer from the fluid to the interface is equal to the heat transfer from
the interface to the solid, or vise versa. Considering the following regions as seen in

20
 2. Theory

Figure 2.7, conservation of energy dictates the heat transfer by the equation[10]:

q̇0 = q̇1 + Su (2.14)

Where,
q̇0 is the heat flux from the fluid to the interface
q̇1 is the heat flux from the interface to the solid
Su is any user specified heat flux

Figure 2.7: The boundary interfaces between a fluid (left) and a solid (right) in
StarCCM+[10].

CHT is widely used in CFD for thermal management problems, climate simulations,
electronics simulations and multi-phase problems. Electric machines have numerous
solid materials along with air and a cooling liquid. The heat losses for different
components are all calculated using electro-magnetic simulations and these varying
losses must be applied to different regions. Also, the introduction of thermal and
contact resistances is only possible through interfaces and hence for this thesis, CHT
CFD is a vital tool.

21
2. Theory

22
 3
Methods

This chapter will discuss the various methodologies and modelling techniques used in
this thesis. The various approaches, why they have been implemented, how it differs
from previous thesis work, and how the results affect this thesis will be elaborated
in the coming sections. The boundary conditions, operating parameters, material
properties and physics models used will also be explained here.

3.1 CAD Modelling

Figure 3.1: CAD model of the ERAD motor, created for this thesis.

As the thesis is carried out both at Chalmers and Volvo Car Corporation, using
sensitive geometries from VCC will be against the Confidentiality Agreement that
has been signed between the participants of this thesis. Hence, a new CAD model
is created for the same machine, for the sake of this thesis. This new CAD model,
as shown in 3.1, is the property of Anthony Jayanath Vivek . The Electrical Power
Engineering Department at Chalmers has purchased the motor in question and have

23
3. Methods

used it for various tests on the rig. This motor is then disassembled and the various
components are placed in Chalmers. These components are used to measure the
required dimensions and recreate in CAD, using CATIA V5.

Apart from this CAD, there is another CAD of the machine created by Anil Kumar,
who had previously conducted a similar thesis[2]. This CAD can also be used in
this thesis but has some slight variations in the geometry. The table below shows
the dimensions of Anil’s CAD vs the actual CAD.

Table 3.1: Differences in dimensions from actual ERAD motor and Anil’s geometry.
ERAD motor Anil’s geometry Relative difference
Component
dimensions (mm) dimensions (mm) to ERAD motor
Length 135 129 -4%
Stator outer diameter 210 212 +1%
Stator inner diameter 145 152 +5%
Rotor outer diameter 143 150 +5%
Rotor inner diameter 48 48 0%

As seen above above, the rotor and the stator volumes vary by 6% and 10% re-
spectively and will affect CHT simulation results. Hence a new model is created in
CATIA V5 with the dimensions measured from a motor which was used in an actual
car. Apart from the dimensions, a few other parts are added in the CAD which are
believed to affect simulation results. They are:

• Different modelling of windings explained in Section 3.10 and 3.11.

• Slot liners are added to the windings and are considered vital as their thermally
insulating properties vary temperature profile through the slots.
• End windings liners and high voltage cables are modelled around the end
windings. They are added as they affect heat transfer from air to windings or
vice-versa.
• Glue layers around the magnets are added as their thermally insulating prop-
erties change temperature gradient between the magnets and the rotor.
• The shaft has separate and named surfaces where bearings and gearbox losses
can be applied.
• The air-gap does not consist of smooth cylinders since it has slots on the stator
side.

Any geometry created for simulations will require surface cleanup for StarCCM+
to be able to work with the geometry. For this, the CAD is exported as .stl or
.stp files and these files are then cleaned in a tool named ANSA [22]. ANSA is a
pre-processor used to clean surface geometry, manipulate geometry, create a surface
mesh and also a volume mesh. By using ANSA, identification to surfaces or regions
can be assigned, and this is useful for automation of large simulations.

24
 3. Methods

3.2 Rotational Analysis

The rotor, as the name states, involves rotation, which must be simulated in CFD.
This rotational simulation can be applied in StarCCM+ in the following ways:

• Moving Wall(MW) Approach

• Moving Reference Frame(MRF) Approach
• Sliding Mesh(SM) Approach

In this thesis, for the rotor, simulations are conducted using all 3 approaches on a
test rotor. In case of a simple cylindrical geometry, a MW approach would be the
fastest and simplest way to simulate the rotation. Since the rotor has holes, the
MW approach will not be valid for the walls of the holes, which are not tangential
to the direction of rotation. Below, the effect of applying a MW approach to the
holes of the test rotor can be observed.

Figure 3.2: Velocity distribution in the rotor holes using a moving wall approach.

From the above image, it is seen that velocity within the rotor holes are not correctly
estimated. The aberration in the velocity profile at the rotor and magnet region in-
terface is due to the fact that the interface does not have a rotation imposed on it.
Hence, an MRF approach is used to model the rotation of the holes in the rotor.
An MRF zone is applied around the regions of the test rotor. For this, a cylinder
is created larger than the rotor and then subtracted from the rotor to produce the
air surrounding the rotor. A relative motion is applied to this region as per the
StarCCM+ manual[10].

The velocity of the rotor around the holes is much more uniform in this case. A
comparison for the velocity scalar for the approach with MW and a similar case with

25
3. Methods

MRF is shown in Figure 3.3. Computational time has increased quite significantly
.It takes 30% more time with the same amount of cores to solve an MRF case for
this rotor or a similar geometry. The results are much better on non-tangential walls
as seen around the rotor holes.

Figure 3.3: Velocity distribution in the rotor holes using a moving wall approach
(left) vs. a similar simulation with an MRF approach (right).

For an MRF approach, the flow estimation between the MRF and the stationary
fluid region is dependant on the calculations at the interface of the MRF and the
stationary region. In this case, for the air-gap between the rotor and the stator, the
thickness is very small (0.8 - 1mm). For the MRF region to cover the rotor com-
pletely, it must envelop the air-gap partially but must not cover the stator region.
Having an MRF region and a stationary region within a gap of 1mm, along with
prism layers to provide good interaction between the two fluid regions, as well as
a finely discretized mesh to capture the Taylor Vortices (which will be discussed in
Section 3.8) would be quite hard in CFD. This would require about 5 layers of prism
layers on the interface and the solid walls, meaning a Low-y + wall treatment, as well
as at least 10 polyhedral cells to capture the vortices. Hence, the mesh will be very
fine and computationally very demanding.

To avoid this, a hybrid MW-MRF approach is implemented to simulate motion

in the rotor. This approach is inspired from [23] where a combination of MW and
MRF is utilized. The MW approach is implemented on walls having surfaces tangen-
tial to the rotational direction. For the rest of the surfaces, an MRF region is created.

26
 3. Methods

3.3 Mesh Sensitivity Analysis

When conducting any computational analysis involving discretization of elements
(Finite Element Methods, Finite Volume Methods), it is very important to conduct
a mesh sensitivity analysis. A mesh independent solution is preferred as the mesh
can be eliminated as a factor for any discrepancies in the solution. The drawback for
certain cases is that to achieve mesh independence, the domain may need to be dis-
cretized very finely. This increases computational time and memory requirements.
Hence, a balance has to be found to obtain reasonable mesh independence without
needing a fine mesh. In this thesis, a mesh sensitivity analysis is conducted for both
the solid and the fluid domain.

For the fluid domain, the k-ω turbulence model is used along with the All y + wall
treatment. Here, the recommended prism layer growth rate (1.3) and 5 prism layer
numbers on the surface are maintained in all cases for the mesh sensitivity analysis.
A polyhedral mesher with automatic surface repair is used in this case. The mesh
cell size is initially set to a very coarse value of 10 mm and gradually decreased all
the way to 1 mm. The number of cells in the fluid domain respectively increases
from 115,000 cells to 1.9 million cells.

Figure 3.4: The increase in mesh quality for the fluid domain from 0.1 million
(left) to 0.7 million (centre) to 1.9 million (right) cells.

Similarly, for the solid domain, a simple conduction case in a cylindrical rod is setup
in CFD and validated against empirical formulations[24]. The cylinder parameters
are shown in Table 3.2.

Table 3.2: Operating parameters for sensitivity analysis.

Parameters Values
Radius r (mm) 50
Length l (mm) 300
Initial rod end temperature T0 (K) 373
Ambient temperature Tamb (K) 300
Ambient HTC h (W/m2 − K) 25
Rod thermal conductivity k (W/m − K) 50

27
3. Methods

On one end of the rod, a temperature of 373 K is applied and an ambient HTC of
25 W/m2 − K is considered throughout the surface of the rod. For the given rod,
the temperature on the other end due to conduction in the rod, and convection to
the environment, can be calculated as:

T0 − Tamb
Tend2 = Tamb + (3.1)
cosh(αl)

Where,
, T0 , l are defined in Table 3.2
Tambq
α = 2h rk

The calculated temperature is 335.72 K. A simple rod with the same dimensions
as above is created in StarCCM+ and the material properties and boundary con-
ditions are imitated in the CFD software. Different mesh types and mesh sizes are
tested to see the temperature variation. A coarse and fine mesh is tested for a poly-
hedral and a trim mesh. The meshes are shown below.

Figure 3.5: The refinements made for the polyhedral mesher (left) and for the trim
mesher (right) for the simple conduction case using a rod.

3.4 Basic EM Model

As this thesis implements a step-by-step approach towards complexity, the first step
in modelling the motor is to start off with a representative basic EM model. This
consists of concentric cylinders with similar dimensions for the different regions as
seen in Figure 3.6.

28
 3. Methods

Figure 3.6: The basic EM geometry with some dimensions.

A fluid domain is placed around the solid cylinders and the domain extends much
larger than the actual casing in the case of the ERAD motor. Since this is a very
early simulation in the thesis, not much consideration is placed in the size of the
fluid domain. The rotor and the magnets are simulated to have rotation by applying
an MW boundary condition. This domain is seen in Figure 3.7.

Figure 3.7: The basic EM geometry with the fluid domain.

This model is used to compare results to a similar Basic Model used in the previous
Master’s thesis on the same topic[2], at similar operating speeds and input torque.
The boundary conditions for this model is as follows:

Table 3.3: Operating parameters and heat loss input conditions for the simple
case.
Losses in W
Description Operating RPM
Rotor Yoke Magnets Active Windings Stator Yoke
Value 4000 20 5 50 25

29
3. Methods

The physics settings activated for the solid and fluid domain are listed in Appendix
C in Table C.1. The solids and air domains are meshed with a polyhedral mesher
with imprint between contacting surfaces to maintain conformal meshing at the
interfaces. This case is then simulated until the residuals reach a reasonable value
of convergence and the results are analyzed.

3.5 Thermal Resistances

In the ERAD motor used in this thesis, there are several components which are
relatively minute as compared to its other components but are still vital to the heat
transfer in the machine. Some examples of such components include:

• The glue used to attach the magnets to the rotor yoke

• The slot liners insulating the active windings from the stator yoke
• The insulating material around the end windings

These materials, along with many other components with simple geometry (high
voltage cables, insulation at the ends of the stator yoke) can be modelled as thermal
resistances in CHT CFD. These thermal resistances replace the intended compo-
nent, meaning one of the surrounding components must be altered in dimensions
to replace the space created by removing this component. By doing so, there is
a possibility of changing the thermal gradients in this component, but considering
that the components that are replaced are very small in dimensions (usually 0.1 to
0.5 mm in this thesis), the effects can be neglected.

Although the dimensions of these components may be very small, since most of
the components considered here are very good thermal insulators, they provide a
considerable thermal resistance to heat transfer. This thermal resistance is calcu-
lated by the Equation 3.2.

L
Rth = (3.2)
k
Where,
L is the thickness of the insulating material in m
k is the thermal conductivity of the material in W/m − K
Rth is the equivalent thermal resistance of the material in m2 − K/W

For the glue holding the magnets to the rotor yoke, the thickness of this glue is 0.25
mm along the length of the magnet and 1.25 mm along the width of the magnets.
Due to such small dimensions, to maintain a conformal mesh across the different
regions, a very fine mesh is required across the glue, the magnets and the rotor at
these regions, as seen in Figure 3.8. Considering the thermal conductivity of the
high temperature epoxy glue to be 0.2 W/m-K, the equivalent thermal resistance is
then 0.00625 m2 − K/W and 0.00125 m2 − K/W across the widths and the lengths
of the magnet respectively. By eliminating the glue from the geometry, the rotor
yoke is increased in volume to replace the glue, and the absence of a thin geometry

30
 3. Methods

allows for a much coarser but still a conformal mesh, as seen in Figure 3.8. With
similar boundary conditions as the first setup, a thermal resistance is applied across
the magnet-rotor interfaces to substitute for the glue. Similarly, cases for thermal
resistances are simulated for the active windings and the end windings, and will be
explained in the coming sections.

Figure 3.8: The difference in mesh quantity without a glue region (left) and with
a glue region (right).

3.6 Contact Resistances

When two solid components are pressed together, depending on the surface finish,
micro hardness, contact pressure and other such material properties of the two
components, a junction is formed with varying contact areas[25] as seen in Figure
3.9.

Figure 3.9: The varying contact areas in a typical junction[25].

The thermal path is majorly conduction through these junctions and only a small
fraction is through convection or radiation in the gaps. The quantity of junctions
with respect to the total surface area of the contacting components is directly propor-
tional to the amount of the heat transferred in the two bodies. Hence the estimation
of contact resistances is very essential in determining the heat transferred in bodies.

31
3. Methods

In an electric machine, there are many components pressed together or assembled

together. Some of these components are:

• Rotor and Shaft

• Shaft and Bearing
• Bearing and Housing
• Stator and Cooling Jacket

The warmest regions of the electric machine are the windings (end and active wind-
ings) or the magnets depending on the operating range and torque of the machine.
The temperatures may exceed 373 K in extreme cases. For the windings, increasing
the temperature would mean decrease in electrical conductivity and hence increase
in copper losses. The Figure A.2 shows the thermal conductivity of copper with
respect to temperature.

For the magnets, they start losing their magnetic properties from around 403 K
and get permanently demagnetized above 423 K. Hence, cooling electric machines is
very important to keep temperatures down in the machine. The rotor, which houses
the magnets, is usually cooled by conduction of the generated heat to the shaft,
which conducts to the casing and to the cooling jacket in turn. The stator is cooled
by directly transferring its heat to the cooling jacket attached to it. The estimation
of the contact resistances between these two components are vital: stator and the
cooling jacket, and the rotor and the shaft. Temperature estimation between the two
interfaces is very crucial to the accuracy of LPN models or even in CFD simulations.

The sensitivity analysis conducted in this thesis includes two thermal contact re-
sistances (TCR) and the gearbox loss in the machine. The TCRs are calculated
based on [26] . The equations and their descriptions are given in Section B.1. The
values of the TCR are shown in Table B.1 and B.2. For the values of resistance,
a sensitivity analysis of the following percentages of the original values are setup:
10%, 30%, 50%, 70%, 100%, 130%, 160% and 200%.

These resistances are added between the respective nodes in the existing LPN of
the ERAD motor. This LPN has been taken from the VeHiCLe project[27]. The
nodes are between the stator and the casing (Resistance R3) and rotor and the shaft
(Resistance R7). There is a gearbox heat loss map also in the LPN model. This loss
has also been added as a parameter in the sensitivity analysis. The gearbox losses
are scaled as a percentage of the original map as: 70%, 80%, 90%, 100%, 110%,
120%, 130%.

These three parameters (R3, R7, and Gearbox Losses) are then setup as a Design of
Experiment (DOE) optimization case in GT-Suite in the LPN model. A total of 448
cases are run for 8 resistances each for R3 and R7 and 7 cases for the gearbox loss
maps. These cases are run for the following operating conditions for the ERAD mo-
tor shown in table 3.4. Test data for these conditions are already available from the
rig measurements conducted for the ERAD motor at Chalmers. The EW temper-

32
 3. Methods

atures recorded from the rig measurements were an average of 329.6 K over 3 sensors.

Table 3.4: Operating Parameters for sensitivity analysis.

Parameters Values
Speed (rpm) 6279
Torque (N-m) 27.9
Coolant Inlet Rate (L/min) 3.8
Coolant Inlet Temperature (K) 301.5

3.7 Material Properties as a Function of Temper-

ature
The physical properties of any material is dependant on the temperature of the
material. When it comes to solids, on varying temperature, density does not vary
greatly in comparison to its base density. This is because solids have very high
densities. Hence, for solids, thermal conductivity is the main property that needs
to be considered as a material property which varies with respect with temperature.
The materials that are considered to have properties varying with temperature are:

• Shaft steel
• Rotor and stator lamination steel
• Winding copper
• Aluminium in cooling jacket

Similarly, for the two fluids in the machine, air and coolant, many properties vary
with respect to temperature. There are four main properties which affect CHT CFD
simulations. These are the density and dynamic viscosity, which are vital to fluid
motion calculations, and thermal conductivity and specific heat, which are mainly
used in thermal calculations. The Section A.2 In Appendix A show the variation of
thermal conductivity of the different materials with respect to temperature.

For each material property, a polynomial fitting curve equation can be extracted
from MS-Excel for the data points. This fitting curve equation is then entered in
StarCCM+ as the input data for the respective material property.

3.8 Air-Gap Modelling

The air-gap in an electric machine is the region between the rotor and the stator
and is a vital geometrical constraint in the machine. The gap exists to allow free
rotation of the rotor without contact with the stator. The larger the gap, the easier
it is to manufacture the motor. The rotor and stator do not need to be manufac-
tured to high tolerances and they will not require high precision balancing. On the

33
3. Methods

other hand, for the electromagnetic inductance through the materials to require the
minimal amount of electrical energy to create the magnetic field, it requires as low a
magnetic reluctance as possible [11]. Since the air has a large magnetic reluctance,
the smaller the air-gap, the more efficient the machine is. Hence, it is a balance
between manufacturability and efficiency for the machine.

Accurate modelling of thermal properties of an EM is very critical to creating good

1D models for vehicle performance characteristics estimation. Solids estimation is
relatively easy as it is purely conduction. Approximation of surface HTCs for solid-
fluid interfaces are the critical parameters [28] as they can be determined by CFD
or experimental results mainly. Creating empirical formulations from these results
can reduce resources required for CFD or experimental setups. An experimental
setup was created by Becker and Kaye to formulate the HTC as a function of the
rotor geometry and fluid properties as per the Equations B.2,B.4,B.5 and B.6 in the
literature [29]. The study showed four different types of flow based on the modified
Taylor’s number T am [28]:

1. Laminar flow for T am < 1700

2. Laminar flow plus Taylor vortices for 1700 < T am < 104
3. Turbulent flow plus Taylor vortices for 104 < T am

Based on the flow type, the Taylor’s number can be estimated as per Equation B.2,
and then the Nusselt number is calculated as per Equations B.4, B.5 and B.6. Fi-
nally, the HTC is calculated based on Equation B.3. These equations are present
in Appendix B in Section B.2. From literature, there is an existence of cross-flow
vortices (swirls in the flow) in the air-gap annulus. These vortices cause an uneven
temperature distribution across the annulus as seen in Figure 3.10.

Figure 3.10: Taylor vortices as seen in the axial cross-section of an EM

(left)[28] and uneven temperature distribution over the stator surface due to Taylor
vortices[30].

34
 3. Methods

Figure 3.11: The mesh in the air-gap with the rotor and stator also shown.

To capture the Taylor vortices that are spoken about above, simple cylinders are used
in StarCCM+ and a very fine mesh around the cylinders is created in the air-gap
annulus. For this, a volumetric refinement zone is incorporated only in the annulus.
The number of prism layers required can be estimated based on the recommended
values for the turbulence model k-ω with an All-y + wall treatment. The prism layer
height is calculated as per Equation 2.12. The mesh in the air-gap zone is shown in
Figure 3.11. A representative stator and rotor with an air-gap width of 0.5 mm is
used. The annulus has 7 prism layers on either side with a growth rate of 1.3. This
case is run for different RPMs corresponding to the different regions of the Taylor
flow.

3.9 Solids Modelling

For the solids in the machine, namely the laminations, shaft, windings, magnets,
insulation and cooling jacket, the CFD software uses the conjugate heat transfer
(CHT) model to estimate heat transfer between solids and liquids. This model can
be used for pure conduction as well. The solids modelling is explained in Section
3.3 under solids mesh sensitivity analysis.

The contact resistances between two solids and thermal resistances are spoken in
detail in the Sections 3.6. As for the interfaces, an Imprint operation is created
in StarCCM+ and this operation creates an intersected surface on two overlapping
surfaces[10]. The surfaces are then remeshed accordingly so as to preserve the geom-
etry of the new surfaces. The mesh on both sides can be conformal or non-conformal.
A conformal mesh produces same number of nodes on both sides of the interface.
This ensures proper mapping of both surfaces, but the downside is that when a large
body is in contact with a complex shape, or a small geometry, the larger body gets
remeshed with a fine mesh similar to the other side, thereby increasing the number
of cells. In contrast, a non-conformal mesh allows unequal number of nodes on both
sides of the interface. This way, the numerical data gets proportionately migrated
to the neighbouring nodes based on the contact ratio. This allows for a coarser
mesh on one side and a finer mesh on the other side, or even to have different mesh

35
3. Methods

types. The disadvantage however is that sometimes, the interface is not completely
mapped on both sides, leading to improper transfer of data across the interface[10].

The different solids in the EM are: shaft, electromagnetic lamination steel (ro-
tor and stator), permanent magnets, copper windings, electrically insulating epoxy,
slot liner, high voltage cables and coolant jacket. The supplier for the ERAD motor
has not provided detailed material properties, and has only provided the material
name or composition. Hence, most properties are obtained from the internet based
on the composition, or from similar supplier data available for generic motors. The
properties for the different materials are listed in Table B.3 in Appendix B in Section
B.3.

3.10 Active Winding Modelling

The active windings in an electric machine, are the copper cables running through
the slots in the stator. They create the electromagnetic field when a current is passed
through them and cause motion in the rotor in an electric machine. At high loads,
they produce a lot of thermal losses due to Ohmic heating, and are also known as
Joules losses as calculated by the Equation 2.1.

Here it is observed that in a winding slot, the amount of copper, or the cross-
sectional area of copper is very crucial to determining the resistance in the cables.
The slots are also filled with epoxy and a slot liner encapsulates the windings from
the stator laminations. Epoxy and the slot liner are very bad thermal conductors
and hence prove to be very vital in thermal calculations[31]. The volume fraction of
copper and epoxy in the slot is close to 0.45:0.55. Another factor to be considered
in this thesis is the contact area of both copper and epoxy with the slot liner. Since
windings are seldom similar in the slot, one can only stochastically or approximately
estimate this value.

3.10.1 Windings Models

Figure 3.12: Windings model from the baseline model[2] and actual windings from
the ERAD machine.

36
 3. Methods

In the previous thesis conducted by [2], the windings were encapsulated in a volume
of epoxy, and this epoxy was the only material in contact with the stator. No slot
liner or thermal resistance was incorporated[2]. Figure 3.12 shows a comparison of
model in [2], which will henceforth be referred to as the baseline model, along with
an image of the actual windings in the ERAD machine.

In this thesis, new models have been proposed to investigate the effects of having
copper in contact with the slot liner, or to be comparatively similar to a complex
windings model. The models in order in Figure 3.13 are known as Complex, Baseline,
Discretized and Anisotropic model. The complex model is a random distribution
of copper cables in an epoxy filled volume. This model is the closest to reality but
due to its complex geometry, surface preparation, meshing and simulation will be
challenging and time consuming and hence is not considered for the final simulation.
The baseline model is similar to the model in [2], but differs slightly in geometry.
The discretized model tries to imitate the complex model where windings are in con-
tact with the slot liner, but is much more simple in geometry, for ease of handling.
Finally, the anisotropic model consists of an equivalent material in the slot volume,
having the combined properties of copper and epoxy. This model will be explained
in detail in Subsection 3.10.2.

Figure 3.13: The different windings models used in this thesis with the legend on
top. L-R: Complex; Baseline; Discretized; Anisotropic.

3.10.2 Anisotropic Material

In the electric machine, there are a few components which have different materials
in the same region. For example, the windings and the epoxy impregnation run in
parallel throughout the slot. Similarly, the lamination in the rotor and the stator
are pressed together with electrically insulating epoxy. The Figure 3.14 shows the
configuration of different material in the EM. Hence, the epoxy and windings in the
slot be considered to be one equivalent material can having different material prop-

37
3. Methods

erties in different directions. This is basic definition of an anisotropic material[32].

Figure 3.14: The thermal resistance configurations of different components in the

EM.

Thus, this equivalent material can be considered to have thermal resistances of

epoxy and copper in a parallel combination in the axial direction of the motor. The
equivalent thermal resistance is then given by the equation B.7[33]. The equiva-
lent thermal resistance in the radial direction is then considered to be in a series
combination and is given by the equation B.8. Similarly, for the rotor and stator,
the thermal resistances of lamination and epoxy to are regarded be in series in the
axial direction and to be in parallel in the radial direction. Then the equations can
be written as B.9 and B.10. These equations and their descriptions are present in
Appendix B in Section B.4.

3.10.3 Windings CFD Modelling

Looking at the three models that can possibly be used in the final CFD simulation,
all three are modelled in StarCCM+ with the same geometrical dimensions and
boundary conditions. This way, the differences in behaviour of the three models
can be observed. For all three models, only one slot out of the 48 in the stator is
considered. The stator length is reduced from 133 mm to 25 mm. The thermal
distribution along the axial direction in an electric machine is quite uniform[34] and
hence, the length axially can be reduced.

For each case, there is the stator, the slot liner, the epoxy and the copper windings
domain. For the anisotropic model, the epoxy and copper are combined to form an
anisotropic material. The front view of all the models are seen in Figure 3.13. A
similar case is run for the three models in consideration, wherein the slot liner is
replaced with an equivalent thermal resistance given by the Equation 3.2. For a slot

38
 3. Methods

liner thickness of 0.3 mm, having a thermal conductivity of 0.143 W/m − K, the
thermal resistance of the slot liner is calculated to be 2.098 ∗ 10−3 m2 − K/W . The
boundary conditions for CFD are shown in the Figure 3.15.

Figure 3.15: The boundary conditions for windings (left) along with a similar case
with a thermal resistance (right).

As seen in the image, the stator has a symmetry wall condition on the side walls.
The top surface has a HTC of 350 W/m2 − K to imitate the cooling jacket effects.
The bottom surface has a HTC of 90 W/m2 − K, which is an average value of HTC
as seen in the air-gap. The ends of the stator have a HTC of 20 W/m2 −K. The case
where the slot liner is replaced with a thermal resistance, has the same boundary
conditions. Now, for each model, there are two different setup for simulations. In
the first case, a heat input of 40W is specified on the windings. The second case is
where a temperature of 360K is imposed on one end of the windings. The first case
shows the heat transfer in the radial direction for each model while the second case
will help visualize the heat transfer in the axial direction.

3.11 End Winding Modelling

The end-windings (EW) are the hottest region of the electric machine for most op-
erating points[35]. They do not contribute to the electromagnetic to mechanical
energy conversion in the motor, and the thermal losses are due to copper losses.
They are usually surrounded by air, unless they are potted with epoxy or other
insulating material. Since air is a bad thermal conductor, the EW get very warm
during the operation of the EM. Excess heating of the EW would lead to melting
of the epoxy insulation and hence short-circuiting the three (or more) phase high-
voltage (HV) current supply. Hence, cooling of the EW is vital to an EM. The heat
from the EW can be transferred to the air surrounding the EW. If there is sufficient
cooling from the cooling jacket, the stator and active windings would be well cooled

39
3. Methods

and hence the temperature gradient between the active and end windings would lead
to good heat transfer away from the EW to the active windings and to the cooling
jacket.

Similar to the active windings modelling, in this thesis, there are three models
for the EW. The three models are: baseline, anisotropic and discretized model. For
simplification and ease of simulation, only a 450 section of the EW is considered.
Unlike the active windings, the EW simulations are CHT CFD cases, wherein there
is an air domain around the solids, to imitate the actual air in the casing. The EW
are encased in an insulating slot liner like material on both sides of the machine. On
one side, a 3-phase HV cable runs along the complete EW length and the cables have
a very thick rubber casing and would hinder heat transfer drastically during short
operations. On the other hand, for continuous operations of the EM, the HV cables
would be at a very high temperature and would add to the heat in the EW. For
each model, there are two cases, wherein one model has all the geometries included
(the insulating liner and the HV cables), and the other case has the liner and the
cables replaced with an equivalent thermal resistance calculated from Equation 3.2.
The Figure 3.16 shows the domains for CFD for both cases.

Figure 3.16: The boundary conditions for EW (left) along with a similar case with
a thermal resistance (right).

For the air domain, all the outer walls have a symmetry condition. The face closest
to the windings are meant to represent the stator and are kept at 325 K. The
extrusion in the domain represents the air-gap and is an air inlet velocity of 1
m/s. The copper or windings is specified a heat source of 25 W. For the case
with the thermal resistances, only the insulating liner has a thermal resistance of
1.75 ∗ 10−2 m2 − K/W and the liner and the HV cables have a combined thermal
resistance of 4.87 ∗ 10−2 m2 − K/W . For the anisotropic material, the equivalent
parallel conductivity (in the direction of the copper cables) is in θ-direction and the
series conductivity is in the radial direction.

40
 3. Methods

3.12 Cooling Jacket & Coolant Modelling

The cooling jacket is one of the most vital components of the EM. As explained in
Section 3.11, it extracts heat from the machine at the stator outer diameter. The
heat dissipated in the rotor due to the machine operation is either conducted to the
shaft or convected to the air in the air-gap. The heat transferred to the shaft is then
conducted to the casing which conducts to the cooling jacket. The heat convected to
the air-gap is transferred to the end windings or the stator, or the alternate occurs,
where the heat from the EW or the stator is convected to the air domain. The heat
from the windings and the stator is conducted to the cooling jacket again. Hence,
the jacket acts as the major source of heat dissipation away from the EM.

A typical cooling jacket around the stator is a spiral or a helical coil. When this
helical coil is unfolded to be a flat surface, this can then be considered to be a flat
plate with channels flowing through it. A project was conducted by FKFS, Germany
on a cooling plate[36]. The project included measuring temperature and pressure
from probes and sensors in the channels of the plate and modelling the flow in CFD
to compare differences in measured and simulated data. This literature is used as
the reference for the modelling of the physics and mesh of the cooling channels in
this thesis.

Figure 3.17: The cooling plate modelled for this thesis.

The Figure 3.17 is setup with the initial boundary conditions for this case and is
shown in Table 3.5.

41
3. Methods

Table 3.5: Boundary conditions for CFD simulation of literature cooling plate.

Boundary conditions Value

Surface temperature of the bottom of the plate Tsurf (K) 333
Inlet temperature of the coolant Tinlet (K) 313
Coolant inlet flow rate (L/min) 3

The coolant consists of a 50% glycol-water mixture and the cooling jacket is an
Aluminium alloy. The material properties for the cooling jacket is available in Table
B.3, and for the coolant in Appendix A in Section A.2. The coolant is modelled as
a material whose properties is a function of temperature, as specified in Section 3.7.
The physics for the coolant and the jacket are specified in Table C.2 in Section C.2.
The prism layer settings are adjusted as per the flow rate and is according to the
Equation 2.12 in the Section 2.4.3.

The Reynold’s number for the flow for different mass flow rates is calculated by
the Equation 3.3. This number then defines the type of flow in the channels. The
flow can be laminar for low Re values (<2300), transitional (2300<Re<104 ) or tur-
bulent for high Re values (>104 ).
ρ ∗ Vinlet ∗ Lc
Recool = (3.3)
µ

Where,
ρ is the density of the coolant in kg/m3
Vinlet is the inlet velocity in m/s
V olume
Lc is the characteristic length of the channel defined by Area of c/s
in m
µ is the dynamic viscosity of the fluid in kg/m-s
For very high mass flow rates, the flow is turbulent, but for low mass flow rates,
the flow is laminar or transitional. A turbulence model is chosen for the cooling
flow for all flow rates because the chosen turbulence model (k-ω SST) can calculate
flow conditions for laminar and transitional flow conditions equally good. Another
reason for choosing a turbulence model for laminar conditions is that even though
the flow may be laminar in the straight channels, the jacket has 180◦ bends and at
sharp bends, the flow tends to be turbulent [13][14][15].

The results of this simulation, are compared against the results from the litera-
ture. The simulation model which produces reasonably good results along with not
consuming too much resources, is chosen for the thesis. This simulation’s mesh and
physics setup is used for the ERAD motor cooling jacket.

3.13 Complete Model & Validation

Once most major components having dissimilar physics are modelled on a compo-
nent level, their physics and mesh settings are implemented into a full EM model.

42
 3. Methods

The boundary conditions for the full EM simulations require the losses for active
components, the coolant conditions and thermal resistances between certain compo-
nents. The losses are obtained by conducting electromagnetic simulations. This is
done in ANSYS Maxwell, and this thesis is using loss maps obtained from [3]. The
loss maps are provided for different regions of the machine: Rotor, stator, magnets,
copper losses and total losses. The copper losses consist of losses in the active wind-
ings, the end windings which form a ring around the ends and the end windings in
the axial direction which connect the ring end windings and the active windings.
The copper losses are then distributed among these regions based on their volume
fractions, which is acquired from CAD. Table 3.6 shows the loss distribution in these
regions. Figure 3.18 is an ERAD motor loss map showing operating speed in RPM
in the X-axis and the operating torque in N-m in the Y-axis. The contours show
the loss levels.

Figure 3.18: Loss maps of copper losses (left) and rotor losses (right)[3][2].

Figure 3.19: Different regions of the end windings and active windings in this
thesis.

43
3. Methods

Table 3.6: The volume fractions of different winding regions.

Volume fraction
Component
%
Active windings 47.6
EW axial (left side) 4.5
EW axial (right side) 4.5
Ring EW (left side) 21.7
Ring EW (right side) 21.7

The thermal and contact resistances are calculated as per Section 3.5 and Section
3.6. Contact resistances are included between the rotor and the shaft, and the stator
and the cooling jacket. Thermal resistances are used to replace the magnet glue and
slot liner.

At Chalmers, an ERAD motor has been tested on a rig setup. The motor was
connected to a DC machine to measure torque output. The controller can manip-
ulate the operating speed and torque output from the machine. Coolant is passed
through the machine at steady temperatures and flow rates and the temperature is
measured at inlet and outlet. Apart from this, a temperature sensor exists in each
phase of the windings (three in total), at the end windings. The machine is run at a
constant speed and torque input until the temperature readings in all sensors attain
a steady state. This temperature is then recorded. For this thesis, one such case is
taken for comparison and the operating conditions are stated below in Table 3.7.

Table 3.7: Operating parameters for the rig setup comparison simulation.

Operating parameter Value

Operating speed (RPM) 2000
Torque (N-m) 60
Coolant inlet temperature (K) 289.77
Coolant inlet mass flow rate (L/min) 2.9
Total copper losses (W) 571.6
Stator loss (W) 104.54
Rotor loss (W) 1.63
Magnets loss (W) 0.16

Apart from the above stated losses, the shaft has bearings, which has losses based
on the operating load on the bearing. This bearing loss is estimated as a function of
operating speed and motor torque as per [3] and [37], which is based on deep-groove
ball bearings. The loss Ploss,Be is given by the equation 3.4 and this loss is applied
as a flux, based on the area of cross-section of the bearing surface on the shaft.
rB e ∗ µB e ∗ T ∗ ωm
Ploss,Be = (3.4)
rRotor_outer

44
 3. Methods

Where,
RB e is the bearing bore radius
µB e is the friction coefficient for the bearing which is assumed to be 0.0015
ωm is the operating speed of the motor in rad/s
rRotor_outer is the outer radius of the rotor

All materials specified in Section 3.7 in this machine are modelled with material
properties as a function of temperature. The anisotropic material is chosen for
modelling the all the windings, as explained in Section 3.10 and 3.11. The physics
for the solids and fluids are as given in Table C.2 and C.1. The mesh settings are
similar to the component level mesh settings and the figure below show the volume
mesh in the full machine simulation.

Figure 3.20: Full ERAD machine volume mesh in the XZ section plane (left) and
the YZ section plane (right).

The simulation is run on the VCC cluster with allocations up to 10000 cores. The
simulation is run with a low Courant number of 0.5 as the mesh quality is reasonably
low enough to provide reliable results, but will not consume much resource for the
simulation. The full simulation setup has about 63 million cells. The cell distribu-
tion is mainly in the air domain in the air-gap and is explained in Table 3.8.

Table 3.8: Volume cell count in different regions of the simulation.

Number of Cells
Region
(x106 )
Air 51.5
Coolant 3.7
Solids 8.2
Total 63.4

45
3. Methods

One simulation from [2] is taken as a case in this thesis to compare the differences in
results. The operating conditions are as given below in Table 3.9. Simulation from
[2] will henceforth be referred to as the Baseline simulation. The main differences
in the two models are as stated below:

• Material properties in the baseline simulation have constant values.

• The geometry varies as mentioned in Section 3.1.
• Contact resistances are not implemented in the baseline model.
• Bearings losses are not implemented in the baseline model.
• The coolant is modelled as a laminar flow in the baseline model.
• The mesh is much coarser in the baseline model.
Table 3.9: Operating parameters from the baseline simulation [2].

Operating parameter Value

Operating speed (RPM) 1000
Torque (N-m) 100
Coolant inlet temperature (K) 313.15
Coolant inlet flow rate (L/min) 3.9
Total copper losses (W) 1227
Stator loss (W) 69.77
Rotor loss (W) 1.21
Magnets loss (W) 0.001

3.14 1D LPN Model

The existing LPN model for an electric machine, taken from [3], is a 9-node model.
The LPN in can be seen in Figure 3.21.

Figure 3.21: The 9-node LPN network from [3].

This LPN considers the whole stator and rotor as one solid lump each. In [3], the
Biot number has not been considered for different lumps. The importance of Bi

46
 3. Methods

has been explained in 2.2 and is calculated as per Equation 2.2. For an LPN to
produce reasonable results, Bi ≤ 0.1 [7]. Approximate Biot number calculations for
the ERAD motor geometry are shown in Table 3.10.

Table 3.10: Calculated Biot number for the ERAD motor.

Geometry Biot number

Rotor at low speeds (<2000 RPM) 0.07
Rotor at high speeds (>2000 RPM) 0.14
Stator 0.31

It can be observed from the table that the rotor is slightly above the acceptable
limit of 0.1 at high speeds and the stator, for all operating conditions, is well above
the limit. Hence, the LPN must be further discretized to have the Biot number of
each lump well below 0.1. To avoid this situation, in this thesis, an LPN model
is proposed to be created based on thermal simulations on simple 2D geometries.
By running 2D CHT simulations in StarCCM+ for a 1/8th section of a rotor for
different conditions, the heat flux directions and the regions of large temperature
gradients in the rotor can be estimated. These regions can then be created as a
lump in the new LPN. Similar boundary conditions can be imitated and the results
are compared. This LPN model is created in GT-Suite in the Integrated Simulation
Environment[38]. In the 2D simulation, there is a 1/8th section of a rotor consisting
of the shaft, rotor, two magnets and one air-pocket. There is a defined temperature
on the shaft surface of 360 K and a HTC at the air-gap surface of the rotor of
90W/m2 − K. The walls on either side of the rotor has a symmetry condition. The
temperature distribution and heat flux is shown in Figure 3.22.

Figure 3.22: A scalar plot of the temperature distribution (left) and a vector plot
of the heat flux (right) obtained from the 2D rotor.

It is seen from the image above that the rotor can be discretized based on major
temperature gradients, or based on heat flux distributions in the simulation. A con-
cept model is created and the simulations are run with similar boundary conditions.

47
3. Methods

The initial concept shows similar temperature distribution in many small lumps,
which implies that those small nodes can be lumped together into much larger ones.
The concept model started as 14-node model, and the next iteration model consists
of a 11-node model, as seen in Figure 3.23.

Figure 3.23: The 11-node LPN model created in GT-Suite for a 1/8th section of
the rotor.

3.15 2D Thermal Model

For the thesis, a 2D thermal model is created for a 1/16th section of the rotor with
a magnet and the shaft. This section includes one magnet, half the air pocket and
a section of the rotor and the shaft. A symmetric condition on the end walls mimic
the whole rotor. The calculations are carried out in MATLAB and are carried out
with a forward difference method (FDM) by an explicit method which is central
in time and forward in space. Basic geometry parameters are input from the user,
defining dimensions and locations of the shaft, rotor, magnet and the air-pocket, in
polar (r −θ) coordinate system. Similarly, the grid is defined in the polar coordinate
system and the calculations are also conducted in the same system. For visual pur-
poses, the geometry and the solution is converted to a Cartesian (XY) coordinate
system.

For a basic setup, a temperature is imposed on one or more boundaries and an

initial temperature is set to all nodes, considered to be an ambient temperature.
In the code, the total time and the time step for the simulation can be defined.
The temperature from the initial time step and the boundary conditions are used

48
 3. Methods

to calculate the temperatures of nodes for the next time step. This is done by using
the neighbouring node values for the calculation of the current node. Hence the ap-
proach is called forward in time and central in space. The final results can display
both the final steady state temperature values in the components, as well as the
temperature distribution in each time step. This helps visualize the heat transfer
paths in the components and can be utilized for designing cooling components.

Figure 3.24: The grid for a 1/16th section of the rotor with the different regions.

49
3. Methods

50
 4
Results

In the chapter, the results of various simulations conducted in this theses is seen,
along with discussions about the results or the methodologies used.

4.1 Mesh Sensitivity Analysis

For the simulations set up in Section 3.3, the results of the mesh sensitivity analysis
is seen for the fluid domain in Figure 4.1. The graph shows the number of cells in
the fluid domain vs. the temperature on the stator’s air-gap surface. It is observed
that after 1.5 million cells, the temperature change is less than 1 degree. Hence the
mesh settings for the simulation with 1.5 million cells in the air domain can be used
for the full machine simulations.

Figure 4.1: The mesh sensitivity analysis for the fluids domain for the number of
cells in the domain vs. the temperature.

For the various solid mesh modelling techniques, the results are seen in Figure 4.2.
For all meshes, the temperature distribution in the rod is very similar. Since a
polyhedral mesh is used for the fluid domain, it is recommended by StarCCM+ user
manual to maintain a conformal mesh in CHT simulations. Hence, the polyhedral
mesh is selected for solids and a mesh size of 3.5mm with refinements at the surfaces
is chosen as the appropriate mesh setting for solids.

51
4. Results

Figure 4.2: The temperature distribution over the rod for different mesh settings
and types.

4.2 Basic EM Model

For the basic EM model, Figure 4.3 shows the Specified y + HTC (for y + = 30) for
the baseline model[2] and the model in this thesis for similar operating parameters
and boundary conditions. Figure 4.4 shows the temperature distribution for the
same two models. From Figure 4.4, it is observed firstly that the geometry differs
greatly for both models. In the basic model, the fluid domain is much larger, and
the model does not have a distinct magnet or end windings region. The relative
difference of results from the basic EM model to the baseline model[2] are seen in
Table 4.1.

Figure 4.3: Specified y + HTC for the baseline model (left)[2] and the basic EM
model (right).

52
 4. Results

Figure 4.4: Temperature distribution for the baseline model (left)[2] and the basic
EM model (right).

Table 4.1: Differences in results for temperature and Specified y+ HTC of the
basic EM model from the baseline model.

∆ Temperature (%) ∆ HTC (%)

Rotor Stator Stator Rotor Stator
outer outer ends air-gap ends
-1.7% +1.5% -12.4% -17.4% +191%

The differences in Table 4.1 highlights the difference in the two simulations. The
air-gap HTC is relatively close to the baseline model and varies mainly due to the
larger fluid domain in the basic EM model.

4.3 Thermal Resistances

For the simulation of the rotor with and without the glue around the permanent
magnets, from Figure 4.5 it is observed that the temperature distribution in the
bodies are similar in both cases. The temperatures on the magnet faces, the glue
faces and the rotor faces are recorded for the first simulation. Similarly, the tem-
peratures across the interfaces are captured and the two results are compared. This
comparison is seen as a difference of the results with the complete geometry against
the one with the thermal resistances, in Table 4.2.

Table 4.2: Differences in results for temperature across the two simulations.

∆Temperature (%)
Magnet Rotor to magnet Rotor Rotor
outer surface upper surface lower surface
0.26% 0.33% 0.26% 0.28%

53
4. Results

Figure 4.5: Temperature distribution in the rotor with glue around the magnets
(left) and with a thermal resistance (right).

The cell count has been reduced by approximately 30 times and in a full machine
simulation, this will reduce computational time drastically. This is an important
reason to implement the listed thermal resistances in the full machine. Apart from
reducing complexity, resources used and time for simulation is also decreased.

4.4 Contact Resistances

The temperature in the end windings is extracted from the LPN model after the
full DOE case is run. These temperatures are then plotted in a 3D map for the
2 resistances and temperature with all 7 gearbox loss cases overlaid in the same
surface plot, as show in Figure 4.6. The numbers of the contact resistances (1-8)
represent the respective values of the base resistance: 10%, 30%, 50%, 70%, 100%,
130%, 160% and 200%.

As observed in Figure 4.6, the change in gearbox losses affect the results quite
drastically. The temperature readings in the end windings are varying by ±5◦ C.
Conversely, changing the resistances from 10% to 200% sees only a minor variation
in the temperature.

54
 4. Results

Figure 4.6: Results of the DOE with varying R3 and R7 values and different
gearbox losses inputs.

Considering only the original gearbox case (gearbox loss map of 100% value) which
is highlighted in red in Figure 4.6 and portrayed in Figure 4.7, for changes in R3,
the temperature gradient is very low. In contrast, by changing R7, the temperature
gradient is slightly more noticeable. This shows that the end winding temperatures
in the LPN model are more sensitive to the resistances between the rotor and the
shaft. Again, the numbers of the contact resistances (1-8) represent the respective
values of the base resistance: 10%, 30%, 50%, 70%, 100%, 130%, 160% and 200%.

Figure 4.7: Temperature distribution in the end windings for change in R3 and
R7 for a fixed gearbox loss.

Since the LPN is most likely to be calibrated against test data (Hällered test track)

55
4. Results

and rig measurements (Chalmers), varying its existing parameters will only cause
deviation in the results and hence the results cannot be taken for its face value.
Hence, these results can only be used to have a general idea of the direction of
change for varying the different parameters.

4.5 Air-Gap Modelling

For the air-gap modelling, simulations are conducted for a laminar case and a tur-
bulent case. From Figure 3.10 it is observed that the cross-flow Taylor vortices are
produced for a turbulent flow. For a case run at 12000 RPM, similar Taylor vortices
across the air-gap annulus are produced as seen in Figure 4.8.

Figure 4.8: Velocity vector in the air-gap annulus for 12000 RPM rotation rate,
when looking at an axial cross-section of the stator and the rotor.

The Figure 4.9 shows the velocity vector in the air-gap annulus for a laminar case
of 200 RPM. There are no Taylor vortices present.

Figure 4.9: Velocity vector in the air-gap annulus for 200 RPM rotation rate, when
looking at an axial cross-section of the stator and the rotor.

The existing LPN calculates the HTC in the annulus as per literature[29]. In Star-
CCM+, the HTC is captured using the Specified y + HTC method at y + = 30.
For the two cases stated above, the HTC on the rotor air-gap surface are 8.9
W/m2 − K(200 RPM) and 66.2 W/m2 − K (12000 RPM). The empirical formula-
tions as per Equations B.2,B.4,B.5 and B.6 in the literature [29] calculate the HTC

56
 4. Results

to be 57.8 W/m2 − K(200 RPM) and 139.3 W/m2 − K (12000 RPM). Hence, using
the Specified y + HTC method is not advisable for this case, as the calculation is very
sensitive to the y + value provided and this value needs to be accurately estimated.
Instead, the temperature can be approximated midway over the air-gap annulus as
seen in Figure 4.10 for the 12000 RPM case. It is seen that the temperature is
quite constant over the majority of the air-gap, axially. At the ends, however, due
to entrainment from the annulus, a relatively lower temperature is observed. The
temperature over the surface is averaged and then compared to the temperature
calculated in the LPN at the midpoint of the air-gap. This way, the LPN can be
calibrated or validated against CFD.

Figure 4.10: Temperature distribution over an arbitrary midway surface over the
air-gap annulus for 12000 RPM.

To observe the temperature in the air-gap annulus more closely, for the turbulent
case, a line probe is created along the axial direction of the annulus, and the temper-
ature is monitored on this line probe. In Figure 4.11, the temperature distribution
on the line probe, over the axial length of the air-gap is shown.

Figure 4.11: Temperature distribution over a line probe over the air-gap annulus
for 12000 RPM.

The temperature varies accordingly to the velocity in the annulus at that point, and
hence can be deduced as a result of the Taylor vortices. For this case, with a con-

57
4. Results

stant temperature over the rotor air-gap surfaces, a small gradient in temperature
can be expected, but when there are heat sources over magnets, windings and the
laminations, a larger variation in temperature may occur. This is something which
can be considered for the LPN.

This modelling case is conducted mainly for simple smooth cylinders, but in the
case of the ERAD motor, the windings in the slots do not coincide with stator inner
surface. These slots are also run as a simple case to see the difference in result. The
velocity vector in the air-gap annulus for 12000 RPM is shown in Figure 4.12.

Figure 4.12: Velocity vector in the air-gap annulus for 12000 RPM rotation rate,
when looking at an axial cross-section of the stator and the rotor.

The graph below shows the temperature distribution in the slot for a line probe
across the axial direction. The X-axis varies from Figure 4.11 because the axis of
the new geometry is different and the length of the stator and rotor was increased
to be more representative of the ERAD machine dimensions.

Figure 4.13: Temperature distribution over a line probe over the air-gap annulus
for 12000 RPM.

4.6 Active Windings Modelling

For the windings models with 40 W heat source on the copper region, the temper-
ature distribution in the different models is shown in Figure 4.14. As observed in
the image, in the baseline model, since the windings are encapsulated by the epoxy,
the temperatures in the windings are much higher than in the other cases. The
temperatures are all representative and do not have the melting point accounted for.

Figure 4.15 shows the temperature distribution for the same case but for models
where the slot liner is replaced with an equivalent thermal resistance. These cases
show similar temperature distribution and with fewer cells in the solid domain.

58
 4. Results

Figure 4.14: Temperature distribution in the different models with similar bound-
ary conditions for active windings. L-R: Baseline model; anisotropic model; dis-
cretized model.

Figure 4.15: Temperature distribution in the different models with similar bound-
ary conditions for active windings with a thermal resistance. L-R: Baseline model;
anisotropic model; discretized model.

Figure 4.16: Temperature distribution in the different models with similar bound-
ary conditions for active windings with a specific temperature imposed at one end
of the windings. L-R: Baseline model; anisotropic model; discretized model.

For the final case with a temperature applied on one end of the windings, the axial
transfer of heat in the windings can be seen in Figure 4.16. Similar to the above

59
4. Results

cases, the windings being encapsulated in epoxy are at a higher temperature and
the stator is cooler.

As seen from the results, the baseline model has it’s shortcomings and the other
two models represent reality more closely. For the sake of lesser geometrical com-
plexity and lesser mesh numbers, the anisotropic model is chosen as the suitable
model for active windings in this thesis.

4.7 End Windings Modelling

For the end windings, the temperature distribution is observed in Figure 4.17. The
baseline and the anisotropic models are at similar temperatures while the discretized
model is slightly cooler.

Figure 4.17: Temperature distribution in the different models with similar bound-
ary conditions for end windings. L-R: Baseline model; anisotropic model; discretized
model.

Another case with similar boundary conditions is run where the liner and the HV ca-
bles are replaced with equivalent thermal resistances. The temperature distribution
for the different models in this case is shown in Figure 4.18.

Figure 4.18: Temperature distribution in the different models with similar bound-
ary conditions for end windings with thermal resistances. L-R: Baseline model;
anisotropic model; discretized model.

The temperatures in the second case are relatively higher than the first case. In
these simulations, there is an air domain, which account for the heat transfer out
of the end windings. By removing certain geometry from the fluid domain and

60
 4. Results

replacing them with thermal resistances, the flow in the domain is altered. This can
be visualized in Figure 4.19.

Figure 4.19: Velocity vector in the fluid domain of the baseline model without
(left) and with thermal resistances (right).

The change in velocities indicate varying Reynold’s number around different surfaces
and hence different HTC on those surface. Thus, to accurately estimate heat trans-
fer, the geometry around fluid-solid interface must be represented as realistically as
possible. As in the active windings, the anisotropic windings model provides rea-
sonable good results with minimal complexity in geometry and a small mesh count
and is hence chosen as the modelling method for the end windings as well.

4.8 Cooling Jacket & Coolant Modelling

For the coolant modelling, various cases are run with different mesh settings and
types. The cases that have achieved mesh independence are presented in this section.
The coolant and cooling jacket are modelled as materials having constant properties
and as materials having properties as a function of temperature. The temperatures
at various sensors are calculated from CFD and are used to compare different models.

It is perceived in Figure 4.20 that the CFD temperatures are a few degrees higher
than literature. The CAD of the cooling plate is based approximately on the geome-
try from literature and is suspected to have sharper bends. These sharper bends lead
to more turbulence and hence more heat dissipated into the fluid[13][14][15]. Hence
the CFD simulation temperatures seem higher than the literature values. This can
be seen in the temperature distribution in Figure 4.21 where the regions in the plate
around the bends are cooler than on the straight channels.

61
4. Results

Figure 4.20: The comparison of the temperatures in different sensors for

literature[36] and the two CFD simulation cases.

Figure 4.21: The temperature distribution in cooling plate.

Although the temperature readings in the CFD simulation with constant material
properties seem to match literature values well, as explained above, the coolant will
be at a higher temperature than in literature. In reality, the coolant operation is
considered for a wide range of temperatures, such as cold starts in winters to a
very hot and sunny day. Hence, the coolant is chosen to be modelled with material
properties as a function of temperature for the full machine.

62
 4. Results

4.9 Complete Model & Validation

For the first case in the complete ERAD motor simulation, the CFD results are
compared against rig measurements for similar operating parameters. These results
can be seen in Table 4.3.

Table 4.3: Results from CFD compared against the rig measurements for similar
operating conditions.

Operating End Windings

Parameters Temperature (K)
Speed Torque Coolant Flow Coolant Inlet Rig CFD
(RPM) (N-m) Rate (L/min) Temp (K) Measurement Calculation
2000 60 2.9 289.77 323.65 323.03

Figure 4.22: The temperature distribution in full ERAD motor for the rig oper-
ating conditions.

The figure above shows the temperature distribution in the machine for the above
operating conditions. The highest temperatures are observed in the end windings
region. For confidentiality purposes, the casing and the cooling jacket cannot be
shown in the image.

For the second case, the simulations for the ERAD motor from this thesis are com-
pared to simulations from [2]. The operating conditions are specified in Table 3.9.
The Table 4.4 shows the difference in temperatures in different regions. It is seen
from Table 4.4 that the baseline model[2] has much higher temperatures in the wind-
ings, as compared to the simulation in this thesis. This is expected from the results
based on the different windings simulations. Possibly, the difference in geometries
can also be a reason for difference in temperature in different components.

63
4. Results

Table 4.4: Temperatures in different regions of the ERAD motor compared between
the baseline model[2] and the simulation from this thesis.

Temperature (K)
Region ∆ T (K)
ERAD
Baseline
simulation
End windings 372 391.8 +19.8
Housing 322.4 319.2 -3.2
Magnets 339.1 318.8 -20.3
Rotor 338.9 326.4 -12.5
Shaft 336.5 322.5 -14
Stator 335.3 319.2 -16.1
Windings 361.1 331.3 -29.8

The Table 4.5 below shows the comparison of the Specified y + HTC at y + = 30
and the reference temperatures for both the baseline simulation and the simulation
conducted in this thesis.
Table 4.5: HTC in different regions of the ERAD motor compared between the
baseline model and the simulation from this thesis along with the reference temper-
atures.

Specified y + Reference
2
Region HTC (W/m − K) Temperature (K)
ERAD Baseline
Baseline ERAD
simulation simulation
Casing 11.3 21.6 351.7 318.4
End windings 17.2 21.1 334.4 368.4
Rotor air-gap 39 35.4 340.9 315.8
Coolant 817.2 300.5 313.2 313.2

Figure 4.23: The temperature distribution in full ERAD motor for the second
case.

64
 4. Results

The Figure 4.23 above shows the temperature distribution in the full domain. Again,
the end windings are the hottest region of the machine. The active windings are
relatively warm as well. The stator and the rotor are relatively cooler. The shaft,
due to high loads, has high bearing losses at one end.

4.10 1D LPN Model & 2D Thermal Model

For a simulation run to steady state in GT-Suite for the setup described in Section
3.14, the temperatures in different nodes are obtained as shown in Table 4.6. The
temperature in certain nodes are averaged to represent the temperature in those
regions. These results are quite similar to the CFD results obtained in Section
3.14. Transient simulations can be run similarly in GT-Suite but since transient
simulations are not included in the scope of this thesis, transient CFD simulations
are not conducted to validate similar simulations from GT-Suite.

Table 4.6: Temperature in different regions from a steady state simulation of the
1D LPN model of the 1/8th rotor.

Region Average Temperature (K)

Shaft 360
Rotor bottom 357.2
Air pocket 350.4
Magnets 344.6
Rotor air-gap
337.9
surface

For the 2D model of the 1/16th section of the rotor, again, transient simulations
are run on MATLAB. The results are plotted on a contour of the rotor section to
analyze the temperature distribution over the rotor. In Figure 4.24, the temperature
distribution over the rotor is shown when a specific temperature condition is imposed
on the left boundary.

Figure 4.24: The temperature distribution in the 1/16th section of the rotor from
the 2D model with the colour bar showing a representative temperature scale.

65
4. Results

66
 5
Conclusions & Future Work

This chapter describes the outcomes of this thesis and the implications of the results
from this thesis.

5.1 Component Modelling

This section explains briefly the inferences from the different component level mod-
elling conducted in this thesis. The new CAD model is more representative of the
ERAD motor in terms of dimensions and also reflects some minor components, which
might have been overlooked otherwise. For the rotational simulations, the hybrid
MW and MRF approach works well without sacrificing computational time or ac-
curacy of the results. Since transient simulations are out of scope, the SM approach
is not considered.

The mesh sensitivity analysis provides a good understanding of the dependency

of the solution on the mesh quality and size. By doing this, the appropriate mesh
can be chosen for both the solid and the fluid domain such that they provide a mesh
independent solution while not compromising on computational resources. Similarly
the basic EM model paints an elementary picture of how to set up a CHT CFD sim-
ulation and how to interpret or visualize the results.

Comparing CFD models to experimental setup, rig data or track measurements

will provide a way to verify the CFD physics and mesh modelling techniques. This
can be seen in the cooling jacket and coolant modelling where CFD results can be
compared to literature with a similar experimental setup. By verifying this CFD
model against this literature, similar modelling can be used for full machine simu-
lations.

Moving on to the different components, thermal resistances show how minor ge-
ometry can be replaced by equivalent resistances to still achieve similar results.
While resistances between multiple solids prove to show good results, the same can-
not be said about replacing geometries at a solid-fluid interface. By doing so, the
fluid domain changes in shape and hence the flow around the body differs. It is
advisable to keep solid geometries around fluid bodies, as this will not only provide
more realistic results, but also reduce cells in the fluid domain, which contains the
majority of the cells in the full simulation.

67
5. Conclusions & Future Work

TCR affect the heat transfer paths in the machine and can act as a tuning pa-
rameter. Although the supplier data may not provide the necessary properties to
accurately estimate the TCR between components, an approximation can be made
and DOE experiments can be conducted to see the effect of each resistance on the
complete system.

Most materials used in this thesis have properties varying with varying tempera-
tures. Some materials show deviation of properties of 20%-40% of the standard
values for temperatures varying within the operating limits of this thesis (290K -
360K). For solids, only thermal conductivity is considered as a function of temper-
ature. For the fluids, the dynamic viscosity, thermal conductivity and density are
considered as a function of temperature. The specific heat capacities of all these
materials are not considered as they are required mainly for transient simulations,
which are out of scope for this thesis.

The importance of modelling the air-gap annulus is explained in this thesis. At

high rotation speeds, the Taylor vortices caused in the annulus makes for an un-
even temperature distribution in the annulus axially. Also, the HTC captured by
CFD software rely highly on the input y + values in the annulus or on the reference
temperature of the fluid. This will vary greatly from empirical calculations of HTC
for the air-gap. Hence, by comparing temperature in the air-gap from CFD esti-
mations and the respective LPN nodes, the LPN models can be calibrated or verified.

For the windings, trying to accurately capture the temperature distribution in the
slot with a realistic distribution of windings and epoxy is very time consuming and
complex. Hence, simplified models are created and their behaviour are observed.
The discretized and anisotropic models show similar results for both the active and
the end windings. The baseline model however shows a higher temperature in the
copper region of the windings as it is encapsulated by a thermally insulating mate-
rial. Since the anisotropic material requires lesser mesh cells, this model is chosen
for the active and end windings.

5.2 Complete 3D Modelling

For the complete 3D modelling, the time taken to setup is almost 1 month. This is
due to intensive geometry preparation and constant update of the geometry so as to
be able to mesh in StarCCM+. The simulation takes less than a week from setting
up physics and boundary conditions to feeding the parameters and finally running
the simulations on the cluster cores. However, the results justify the extensive lead
time. For the case where the CFD simulations imitate rig measurement operating
conditions, there is a good correlation of the end windings temperature, which is
the only available sensor data from the rig.

In the second case where the model is compared to baseline simulations from pre-

68
 5. Conclusions & Future Work

vious thesis work, there exists a large deviation in temperature and HTC values
from the baseline simulation. The difference in different simulation models for both
simulations may be the reason for the variation in results.

This modelling approach can be used in the future for new concept validation or for
calibrating or validating existing or new LPN models. As these CFD simulations
are very time consuming, they can then be used mainly for calibration or validation
of 1D or 2D models. Running complex drive cycle analysis or transient simulations
would not be advisable.

5.3 1D & 2D modelling

The 1D concept rotor model in GT-Suite introduced in this thesis does not rely on
the Biot number of geometries and solely relies on simple 2D CFD simulations con-
ducted on a section of the geometry. These 2D simulations take very little time to
setup and run. These 1D models can then be replicated for the full geometry for the
solids. For the fluids domain, the complex 3D simulations can provide the missing
data,or this data can be taken from literature. These 1D models can then be used to
run drive cycle analysis. Since these models are created based on geometries input
by the user, new concepts can be tested with these models to verify the efficiency
and performance. Hence, this model can prove to be very reliable in concept and
development phase.

Similar to the 1D model, the concept 2D MATLAB model spoken about in this
thesis calculates the thermal diffusion in solids, much like the 2D CFD simulations
in StarCCM+. This model can be implemented as a Crank-Nicholson method as
well with implicit solvers to gain better accuracy. This solver can then be made
as blocks in Simulink and like the 1D model in GT-Suite, can be run for transient
cases. If programmed in Python, this code can then be made to run on parallel
cores and the process can be run on supercomputer clusters.

With the 1D model, there is reasonably good results for very low lead times. It pro-
vides node temperatures but does not describe the temperature distribution within
components. The 3D CFD model on the other hand gives precise results but is
very expensive computationally. The heat transfer paths, temperature variations
and flow fields can be visualized to a great extent. The 2D model combines the best
parts of the 1D model and the 3D model to give good results within reasonably fast
times and at the same time, providing a representative temperature distribution in
the components.

5.4 Future Work

This thesis covers the 3D thermal modelling aspects of an electric machine. There
are many facets that can be continued or improved upon in the same line as this
thesis. Some of them are listed below:

69
5. Conclusions & Future Work

• The 1D LPN model can be developed in conjunction with 3D CFD and vali-
dated simultaneously.
• CFD can be coupled with electromagnetic and mechanical codes in compatible
software such as StarCCM+ Electromagnetics, ANSYS Maxwell and Fluent,
or OpenFOAM.
• The overall mesh quality can be improved albeit at the expense of higher mesh
count and higher computational requirements.
• Inversely, the mesh around the air-gap can be simplified and certain geometries
that have been replaced with thermal resistances can be added back to the
model.
• Contact resistances can be looked into in more detail for all components.
• Gearbox losses can be implemented into the simulations.
• The 2D thermal model can be calculated using implicit solvers such as Crank-
Nicholson methods for higher accuracy and better stability.
• The 1D and 2D models can be created in Python so that the calculations can
be parallelized on the computer cores or to be able to utilize cluster cores,
thereby decreasing computation time.

70
 Bibliography

[1] Volvo Car Corporation, ERAD Motor

www.media.volvocars.com/global/engb/media/photos/188278/drive-e-
electric-rear-axle-drive-erad
(visited on 08 August 2019)
[2] Anil Kumar, (2018) Electric motor internal heat convection modelling and anal-
ysis, Master’s thesis in Applied Mechanics, Department of Mechanics and Mar-
itime Sciences, Chalmers University of Technology
[3] Emma Arfa Grunditz, (2016) Design and Assessment of Battery Electric Vehi-
cle Powertrain, with Respect to Performance, Energy Consumption and Electric
Motor Thermal Capability, Thesis for the Degree of Doctor of Philosophy, De-
partment of Energy and Environment, Division of Electric Power Engineering,
Chalmers University of Technology
[4] Y. A. Çengel, (2008) Introduction to Thermodynamics and Heat Transfer,
Mcgraw-Hill
[5] Georgios D. Demetriades et. Al., (2010) A Real-Time Thermal Model of a
Permanent-Magnet Synchronous Motor, IEEE Transactions on Power Electron-
ics, Vol. 25, No. 2, February 2010
[6] Sonja Tidblad Lundmark and Anders Bergqvist and Svetla D. Chakarova-Käck,
(2017) Coupled 3-D thermal and electromagnetic modelling of a liquid-cooled
IPM traction motor, IEE 978-1-5386-1317-7/17
[7] Wikipedia, Lumped Element Model,
https://en.wikipedia.org/wiki/Lumped_element_model
(visited on 28 January 2019)
[8] Adrian Bejan and Allan Kraus, (2003) Heat Transfer Handbook, John Wiley &
Sons, Inc.
[9] Austin Hughes and Bill Drury, Electric Motors and Drives, (2013) Fourth Edi-
tion, Elsevier Ltd.
[10] Simcenter StarCCM+ v1306, The Steve Portal StarCCM+ Documentation,
https://documentation.thesteveportal.plm.automation.siemens.com/starccm
plus_latest_en/index.html?param=WLn8nauthLoc=https://thesteveportal.
plm.automation.siemens.com/AuthoriseRedirect
[11] Zhao-Dong Xu and Ying-Qing Guo and Jun-Tao Zhu and Fei-Hong Xu, (2017)
Intelligent Vibration Control in Civil Engineering Structures, Chapter 5 - De-
sign and Parameters Optimization on Intelligent Control Devices, pages 151-
171, Academic Press, Oxford, 978-0-12-405874-3.

71
Bibliography

[12] Lars Davidson, (2018) Fluid mechanics, turbulent flow and turbulence mod-
elling, Division of Fluid Dynamics, Department of Applied Mechanics,
Chalmers University of Technology
[13] C. Carlsson and E. Alenius and L. Fuchs, (2015) Swirl switching in turbulent
flow through 90deg pipebends, P hysicsof F luids27, 085112
[14] Leo H. O. and Metodi B. Zlatinov and Alexander J. Smith
and Guongjun Cao, (2013) Turbulent pipe flow downstream of a
90deg bend, Journalof F luidM echanics735
[15] Athanasia Kalpakli, (2012) Experimental study of turbulent flows through pipe
bends, Master’s thesis, Technical Reports from Royal Institute of Technology
KTH Mechanics
[16] Stephen B. Pope, (2012) Turbulent Flows, Cambridge University Press, ISBN:
9780511840531, June 2012
[17] J.E. Bardina and P.G. Huang and T. J. Coakley, (1997) Turbulence Modeling
Validation, Testing, and Development, NASA Technical Memorandum 110446
[18] David C. Wilcox, (1998) Turbulence Modeling for CFD, Second edition. Ana-
heim: DCW Industries, 1998. pp. 174
[19] CFD-Online, K-epsilon models,
www.cfd-online.com/Wiki/K-epsilon_models
(visited on 16 February 2019)
[20] Frank P. Incropera and David P. Dewitt, (2011) Fundamentals of Heat and
Mass Transfer, 7th Edition, John WIley Sons Publication
[21] The Engineering ToolBox, Emissivity Coefficients Materials
https://www.engineeringtoolbox.com/emissivity-coefficients-d_447.html
(visited on 02 February 2019)
[22] BETA CAE Systems, ANSA pre-processor, https://www.beta-
cae.com/ansa.htm
[23] Teddy Hobeika and Simone Sebben, (2018) CFD investigation on wheel rotation
modelling, Journal of Wind Engineering Industrial Aerodynamics 174 (2018),
pages 241–251, Elsevier Publications
[24] Physics Stack Exchange, (2017) Determining the temperature at the end of a
rod when one end is heated
https://physics.stackexchange.com/questions/296858/determining-the-
temperature-at-the-end-of-a-rod-when-one-end-is-heated (visited on 03 March
2019)
[25] Thermopedia, (2011) Thermal Contact Resistance
http://www.thermopedia.com/content/1188/ (visited on 16 May 2019)
[26] K. N. Babu, (2015) Thermal Contact Resistance: Experiments and Simulation,
Master’s thesis in Applied Mechanics, Department of Mechanics and Maritime
Sciences, Chalmers University of Technology
[27] Virtual HybrId CooLing (VeHICLe) https://www.chalmers.se/en/projects/Pages/Virtual-
HybrId-CooLing-QVeHICLeQ.aspx (visited on 11 Aug 2019)
[28] David A. Howey and Peter R. N. Childs and Andrew S. Holmes, (2010) Air-
gap convection in rotating electrical machines, IEEE Transactions on Industrial
Electronics ( Volume: 59 , Issue: 3 , March 2012 ), pages 1367 - 1375

72
 Bibliography

[29] K. M. Becker and Joseph Kaye, (1962) Measurements of Diabatic Flow in an

Annulus with an Inner Rotating Cylinder, ASME
[30] Md. Lokman Hossain and Rebei Bel Fdhila, (2017) Air-Gap Heat Transfer in
Rotating Electrical Machines: A Parametric Study, 9th International Confer-
ence on Applied Energy, ICAE2017, 21-24 August 2017, Cardiff, UK, ScienceDi-
rect
[31] Mariia Polikarpova and Pia Lindh and Chris Gerada and Marko Rilla and
Ville Naumanen and Juha Pyrhonen, (2015) Thermal effects of stator potting
in an axial-flux permanent magnet synchronous generator, Applied Thermal
Engineering 75 (2015) 421e429
[32] Wikipedia, Anisotropy, https://en.wikipedia.org/wiki/Anisotropy (visited on
21 March 2019)
[33] Massachusetts Institute of Technology, Thermal Resistance Circuits, Mas-
sachusetts Institute of Technology (visited on 22 March 2019)
[34] JinXin. Fan and ChengNing. Zhang and ZhiFu. Wang and E. G. Strangas,
(2010) Thermal Analysis of Water Cooled Surface Mount Permanent Magnet
Electric Motor for Electric Vehicle, IEEE, 2010 International Conference on
Electrical Machines and Systems
[35] R. Aziz and G.J. Atkinson and S. Salimin, (2017) Thermal Modelling for Per-
manent Magnet Synchronous Machine (PMSM), International Journal of Power
Electronics and Drive System (IJPEDS), Vol. 8, No. 4, December 2017, pp.
1903 1912
[36] Karras N. and Kuthada T. and Wiedemann J., (2014) "An Approach for
Water Jacket Flow Simulations", SAE Technical Paper 2014-01-0659, 2014,
doi:10.4271/2014-01-0659
[37] SKF, Rolling bearings, https://www.skf.com/binary/138-121486/SKF-rolling-
bearings-catalogue.pdf (visited on 13 June, 2019)
[38] Gamma Technologies, GT-SUITE ISE, https://www.gtisoft.com/gt-suite/gt-
suite-overview/

73
Bibliography

74
 A
Appendix A - Additional Images

A.1 Introduction

Figure A.1: The initial planning Gantt chart.

I
A. Appendix A - Additional Images

A.2 Methods
A.2.1 Material Properties as a Function of Temperature
Solids:
Copper considered to be a normal copper wire used for electrical appliances.

Figure A.2: The thermal conductivity of copper varying with temperature.

Lamination steel properties obtained from:

https://www.tatasteeleurope.com/static_files/Downloads/Construction/
Energy%20and%20Power/Electrical%20steel/non-oriented%20electrical%20steels.pdf

Figure A.3: The thermal conductivity of lamination steel as a function of temper-

ature.

II
 A. Appendix A - Additional Images

Shaft steel is assumed to be AISI1046 steel, which is a commonly used steel for
shafts.

Figure A.4: The thermal conductivity of the shaft steel varying with temperature.

The cooling jacket is assumed to be Al7075-T6 Aluminium alloy which is a sturdy

and lightweight, commonly used alloy for cooling cases.

The slot liner is assumed to be DuPont Teijin Films Mylar® A Polyester Film,
48 Gauge, whose properties are available at:
http://www.matweb.com/search/datasheet.aspx?matguid=3e9bc72bfe08408aa7
acbd568d15a895ckck=1

The permanent magnets are assumed to be Neodymium Iron Boron magnets, which
are very commonly used magnets and their properties are available from:
https://www.eclipsemagnetics.com/media/wysiwyg/datasheets/magnet_materials
_and_assemblies/ndfeb_neodymium_iron_boron-standard_ndfeb_range
_datasheet_rev1.pdf

III
A. Appendix A - Additional Images

Coolant (50% Ethanol - 50% Water):

Figure A.5: The variation of dynamic viscosity with temperature of the coolant.

Figure A.6: The variation of thermal conductivity with temperature of the coolant.

IV
 A. Appendix A - Additional Images

Figure A.7: The variation of density with temperature of the coolant.

Air:

Figure A.8: The variation of dynamic viscosity with temperature of the air.

Figure A.9: The variation of density with temperature of the air.

V
A. Appendix A - Additional Images

VI
 B
Appendix B - Additional
Information

B.1 Contact Resistances

Contact thermal resistance,
1
Rctr = eq ∗θRM S 0.94 (B.1)
( 1.13∗k
σRM S ∗10−6
)( HPc )

Where,
Contact pressure, P MPa
Micro-hardness, Hc MPa
2k1 k2 W
Equivalent thermal conductivity, keq = k1 +k2 m−K
q
RMS surface hardness, σRM S = σ12 + σ22 µm
q
RMS asperity slope, θRM S = θ12 + θ22

The assumed material properties for the various components are listed below:

Table B.1: TCR calculation assumptions and results for rotor and shaft

Component Rotor Shaft

W
Thermal Conductivity m−K 40 51
Roughness µm 10 0.05
Asperity slope 0.25 0.1
Micro-hardness MPa 3800 3800
Contact Pressure MPa 5
keq 44.83
σRM S 10
θRM S 0.27
W
Thermal contact conductance m−K 2672.37
Thermal contact resistance m−K
W
3.74 ∗ 10−4

Similarly for the stator to cooling jacket, a similar TCR is calculated. An air-gap
of 1 micro meter is assumed around this interface as well and the final contact re-
sistance is a sum of the TCR and the thermal resistance due to the air.

VII
B. Appendix B - Additional Information

Table B.2: TCR calculation assumptions and results for stator and cooling jacket

Component Stator Cooling Jacket

W
Thermal Conductivity m−K 40 201.1
Roughness µm 10 0.125
Asperity slope 0.25 0.03
Micro-hardness MPa 3800 1400
Contact Pressure MPa 5
keq 66.72
σRM S 10
θRM S 0.25
W
Thermal contact conductance m−K 3718.94
Thermal contact resistance m−K
W
2.69 ∗ 10−4
Air thermal resistance m−K
W
3.82 ∗ 10−5
Total thermal resistance m−KW
3 ∗ 10−4

B.2 Air-Gap Modelling

0.5
ω ∗ rm ∗ (b − a)1.5
T am = (B.2)
ν

h ∗ Dh
Nu = (B.3)
k
For T am < 1700

2(b − a)/a
Nu = (B.4)
ln[1 + (b − a)/a]

For 1700 < T am < 104

0.367
N u = 0.128T am (B.5)

And for 104 < T am

N u = 0.0.409T a0.241
m (B.6)

Where,
ω is the rotation speed of the rotor in rad/s
b and a are the stator inner and rotor outer diameters respectively in m
ν is the kinematic viscosity of the air-gap fluid in m2 /s
Nu is the Nusselt number
h is the HTC of the fluid in W/m2 − K
Dh is the hydraulic diameter of the annulus in m
k is the thermal conductivity of the fluid in W/m − K

VIII
 B. Appendix B - Additional Information

B.3 Solids Modelling

Table B.3: Material properties for solids in the ERAD machine

Component Density Thermal Conductivity Specific Heat

(Material) kg/m3 W/m − K J/kg − K
Shaft
7800 51 460
(Steel)
Lamination
(Electromagnetic 7800 40 460
steel)
Permanent magnets
(Neodymium Iron 8000 9 400
Boron)
Winding
8900 395 385
(Copper)
Slot
1350 0.2 1700
(Epoxy)
Slot/EW liner
1350 0.143 1700
(BoPET)
HV cable shield
930 0.32 1971
(Insulating rubber)
Cooling jacket
2705 230 900
(7075-T6 Aluminium)

B.4 Anisotropic Material

1 1 1 kC u ∗ AC u + kEpoxy ∗ AEpoxy
= + = (B.7)
Rparallel(Cu,Epoxy) RCu REpoxy L

L L
Rseries(Cu,Epoxy) = RCu + REpoxy = + (B.8)
kC u ∗ AC u kEpoxy ∗ AEpoxy

Where,
L is the length of the copper or epoxy in the slot (both are equal) in m
kC u is the thermal conductivity of copper in W/m − K
AC u is the area of cross-section of copper in m2
kEpoxy is the thermal conductivity of epoxy in W/m − K
AEpoxy is the area of cross-section of epoxy in m2
In this equation, we take in to consideration the volume fraction of copper and epoxy
(from the length and are of cross-section) and hence the equivalent thermal resis-
tance is a function of the volume fraction and the respective thermal conductivity
of copper and epoxy.

IX
B. Appendix B - Additional Information

1 1 1 kSteel ∗ ASteel + kEpoxy ∗ AEpoxy

= + = (B.9)
Rparallel(Steel,Epoxy) RSteel REpoxy L

L L
Rseries(Steel,Epoxy) = RSteel + REpoxy = + (B.10)
kSteel ∗ ASteel kEpoxy ∗ AEpoxy

Where,
L is the length of the copper or epoxy in the slot (both are equal) in m
kSteel is the thermal conductivity of steel in W/m − K
ASteel is the area of cross-section of steel in m2
kEpoxy is the thermal conductivity of epoxy in W/m − K
AEpoxy is the area of cross-section of epoxy in m2

X
 C
Appendix C - Physics Models

C.1 Air and Solids Physics

The table below lists the physics settings used in StarCCM+ for different solids and
the air domain in this thesis.

Table C.1: Physics Continuum for the Basic EM Case

Fluid Continua Solid Continua

Cell Quality Remediation Cell Quality Remediation
Gravity Constant Density
All y+ Wall Treatment Coupled Solid Energy
SST (Menter) K − ω Steady 3D
Coupled Energy
Ideal Gas
Coupled Flow
Steady 3D

C.2 Coolant & Cooling Jacket Modelling

The table below lists the physics settings used for the coolant & cooling jacket
modelling in this thesis.

Table C.2: Physics Continuum for the Coolant & Cooling Jacket Modelling Case

Fluid Continua Solid Continua

Cell Quality Remediation Cell Quality Remediation
Gravity Constant Density
All y+ Wall Treatment Coupled Solid Energy
SST (Menter) K − ω Steady 3D
Coupled Energy
Liquid Polynomial Density
Coupled Flow
Steady 3D

XI