my

.353

Euawfiwm . 3m“?

 

 

 

0. ”W45?“

1..
..

a.
. I.
.

a . V535, 1m. ‘3.an ’bu

 

 

51.335..- 3%

.$¢Sfitt!:1.. 5.5.: .
1.

.l ..
V . k , £
» 1‘: . . II. ..M
t.‘.fil..rtr.nw..unzi JN-g‘vg‘,‘

.. < .

I
(1.. 1. . I I

5.1....

at: .

{a 2..
We. Fudwwv.”
u :3 II 314.:

. .lflikwm.
1... . maxi...
1m .

.fJ‘ “filum
.§wfimfin...rilia
1:1: 2..) 5
2: blast)!!!
(in...

t (3.0:.
x it... .v

3.. : 3.1:!!-

3.75.9...113
5... 3‘0 «
IQ! A

.2...

10:51.2:

It: .3»...-
in

. nTIAE‘
I: 12:21

,h-lt:«l‘ .
.ué!)n)x.t€. 01:33..
1... .A ‘ .

1!. (fill, l: ..n
)6. 13:11.: .I. 3...

 

l “"2“ LIBRARY
237/, Michigan State
University

ca

 

 

 

This is to certify that the
thesis entitled

PISTON DESIGN AND ANALYSIS: PARAMETERIZED AND
COMPLETE FINITE ELEMENT ANALYSIS APPROACH FOR
THE ASSESSMENT OF PISTON PERFORMANCE

presented by

ANDREAS PETROU PANAYI

has been accepted towards fulfillment
of the requirements for the

MS. degree in Mechanical EnflteerinL

 

 

[W40 W

Maflfiofessor’s Signature
6% 25], 2mg

Date

MSU is an Affinnative Action/Equal Opportunity Institution

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 p:/C|RCIDaIeDue.indd-p.1

 

PISTON DESIGN AND ANALYSIS: PARAMETERIZED AND
COMPLETE F INITE ELEMENT ANALYSIS APPROACH FOR THE
ASSESSMENT OF PISTON PERFORMANCE

By

Andreas Petrou Panayi

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Mechanical Engineering

2006

ABSTRACT
PISTON DESIGN AND ANALYSIS: PARAMETERIZED AND COMPLETE FINITE

ELEMENT ANALYSIS APPROACH FOR THE ASSESSMENT OF PISTON
PERFORMANCE

By

Andreas Petrou Panayi

Modeling the thermal and mechanical behavior of a piston is crucial, as it allows
for the prediction of piston characteristics such as dynamics and fiction. These
characteristics directly affect the efficiency of an internal combustion engine. In this
thesis, two studies are conducted; an evaluation of the parameterized piston modeling
method and a development of a finite element program for a full piston model analysis.

For the first study, a parameterized piston model is generated using a set of
selected geometrical parameters from the production piston model according to the
standard of CASE, a cylinder-kit simulation software. A finite element analysis is
conducted for both the parameterized piston model and the production piston model using
the software COSMOSDesignSTAR. A comparison of the two models shows good
agreement in both thermal and mechanical characteristics. CASE is used to demonstrate
the advantages of the parameterized model in the assessment of piston performance.

For the second study, a piston finite element analysis program is developed. The
capabilities of the program are demonstrated using a full piston CAD model. The meshed
geometry of the CAD model is imported and an analysis is performed over a full four-
stroke cycle. The piston’s elastic behavior under thermal and mechanical loads is verified
with COSMOSDesignSTAR and the elastohydrodynamic lubrication analysis results

show that this model is capable of predicting the piston behavior and performance.

Cepyright by
ANDREAS PETROU PANAYI
2006

to my parents...

and to my Darioulla...

iv

ACKNOWLEDGEMENTS

First I would like to give sincere thanks to my advisor, Dr. Harold J. Schock, for
giving me the opportunity to work with him, first as his teaching assistant and then as his
research assistant at the Automotive Research Experiment Station here at Michigan State
University. I would also like to extend my thanks to him for his guidance and faith in me
through my studies and my research.

Special thanks go to Dr. Andre Benard and Dr. Ronald Averill, the members of
my guidance committee. They both provided me with valuable information on finite
element methods, thus strengthening my understanding of the methods and consequently
enabling me to develop the piston finite element analysis program. I would like to thank
Dr. Giles Brereton as well for providing me with invaluable literature references and for
discussing and suggesting solutions to problems I was encountering with numerical
modeling.

I would also like to express my gratitude to Dr. Boon-Keat Chui of Mid Michigan
Research, first for encouraging me to have my first publication and second for all his
guidance in understanding the piston modeling techniques and CASE.

Dr. Mikhail Ejakov of Ford Motor Company deserves special recogtition for
providing me with the piston solid models, as well as for his willingness to discuss the
development of the piston finite element analysis program.

I would like to thank Andrew Fedewa of Mid Michigan Research for sharing his

experiences with CASE and piston modeling and for providing me with the experimental

oil film thickness measurements. Also I thank Mulyanto Poort of Mid Michigan Research
for his invaluable programming tips and advice.

Special recognition goes to the people of the Automotive Research Experiment
Station: Dr. Murad Ismailov, Ed Timm, Tom Stuecken, Josh Bedford, .Iia Ma, Kim Sarbo
and Melissa Jopke for making it a fun workplace.

Also I would like to thank my dear friend Chandra Rome] for sharing our thesis
experiences and racing with me against time.

Finally I would like to thank and express my gratitude to my parents, Petros and
Panayiota, and to my Darioulla, for their love, support and patience, even when I am so

far and away.

vi

PREFACE

Computer simulations are becoming more popular by the day as advances in
electronics decrease computation time. One of the most widely used computational
methods in engineering now is finite element methods. By now this method is well
understood, and because of its robustness it allows for the solution of complex as well as
nonlinear systems. This thesis presents the finite element method as applied to the
analysis of a piston model. Commercial software is utilized, as well as a piston finite
element analysis code developed by the author.

The thesis starts by introducing the basics of finite element methods. The Finite
Element Theory chapter explains the methods used within this work. The aim of this
chapter is to develop an understanding of the theories upon which finite element codes
are built, thus enabling a confidence in using them. The author assumes the reader has
some previous knowledge in calculus, linear algebra, solid mechanics, heat transfer and
fluid mechanics as well as in numerical methods. The theories do not go in depth; thus,
with the assumed knowledge mentioned above, they can be very easily understood.

Also this thesis introduces the use of commercial software, which is a very
important tool for the engineer of today. A parameterized modeling method is presented
which significantly speeds up simulations. A finite element analysis of a parameterized
piston model as well as of a full piston model is conducted with COSMOSDesignSTAR.
A cyclic analysis is performed in CASE.

Furthermore the thesis shows how a custom finite element code can be developed

to meet specific needs. The last chapter clearly explains the steps required in the

vii

development of a piston finite element analysis code. Although the code itself is not
provided, detailed flowcharts are shown.

The images in this thesis are presented in color.

viii

TABLE OF CONTENTS

 

 

 

 

 

 

 

LIST OF TABLES - _ - - - - xi
LIST OF FIGURES - xii
KEY TO SYMBOLS AND ABBREVIATIONS - ..... -- - xviii
CHAPTER 1 INTRODUCTION 1
1.1 Motivation ......................................................................................................... 1
1.2 Previous Efforts in Piston Modeling ................................................................. 2
1.3 Objective ........................................................................................................... 4
CHAPTER 2 THE FINITE ELEMENT METHOD 6
2.1 Introduction ....................................................................................................... 6
2.2 Approximate Solutions ..................................................................................... 8
2.2.1 Galerkin’s Method ............................................................................. 8
2.2.2 Variational Method ............................................................................ 9
2.3 Elements and their Shape Functions ................................................................. 9
2.3.1 Linear Triangular Element ................................................................. 9
2.3.2 Linear Tetrahedral Element ............................................................. 14
2.4 Finite Element Formulation ............................................................................ 21
2.4.1 Heat Equation ................................................................................... 21
2.4.2 Static Analysis ................................................................................. 26
2.4.3 Reynolds Lubrication Equation ....................................................... 37
2.4.4 Incorporation of Boundary Conditions ............................................ 43
2.4.5 Guyan Reduction ............................................................................. 44
2.4.6 Conjugate Gradient Method ............................................................. 46
CHAPTER 3 ENGINE OPERATING PARAMETERS 50
3.1 Introduction ..................................................................................................... 50
3.2 Time Step ........................................................................................................ 51
3.3 Angular Speed ................................................................................................. 51
3.4 Piston Motion .................................................................................................. 51
3.4.1 Sample results .................................................................................. 55
CHAPTER 4 PARAMETERIZATION AND FEA APPROACH FOR THE
ASSESSMENT OF PISTON CHARACTERISTICS 58
4.1 Introduction ..................................................................................................... 58
4.2 Piston model FE analysis.................... ............................................................ 59

ix

4.2.1 Thermal Loading .............................................................................. 61

4.2.2 Thermal Loading Deformation ........................................................ 68
4.2.3 Mechanical Loading ......................................................................... 70
4.2.4 Combined Thermal and Mechanical Loading .................................. 71
4.2.5 Computation Time ........................................................................... 80
4.3 CASE Analysis ............................................................................................... 81
4.3.1 Finite Element Analysis ................................................................... 81
4.3.2 Piston Characteristics ....................................................................... 82
4.4 Experimental Results ...................................................................................... 87

CHAPTER 5 PISTON FINITE ELEMENT ANALYSIS: THE FULL MODEL

 

 

APPROACH - -- 89

5.1 Introduction ..................................................................................................... 89
5.2 Preprocessing .................................................................................................. 89
5.2.1 CAD Model ...................................................................................... 90

5.2.2 Mesh and boundary conditions ........................................................ 91

5.3 Piston Finite Element Analysis (rt-fea) program ............................................ 92
5.4 1t-fea Verification .......................................................................................... 103
CHAPTER 6 RESULTS A - -_ - - -_ _ - ........ - 105
6.1 Introduction ................................................................................................... 105
6.2 Engine geometry and operating conditions ................................................... 105
6.3 Piston geometry and properties ..................................................................... 106
6.4 Mesh Refinement .......................................................................................... 108
6.5 Periodicity and convergence ......................................................................... 112
6.6 Simulation Results ........................................................................................ 114
6.6.1 Temperature ................................................................................... 115

6.6.2 Deformation ................................................................................... 115

6.6.3 Strain .............................................................................................. 116

6.6.4 Stress .............................................................................................. 117

6.6.5 Lubrication Pressure ...................................................................... 118

6.6.6 Oil Film Thickness ......................................................................... 122

6.6.7 Forces and Moments ...................................................................... 123

CHAPTER 7 CONCLUSIONS AND RECOMMENDATIONS - -- - 128

 

 

 

 

7.1 Conclusions ................................................................................................... 128
7.2 Recommendations ......................................................................................... 132
APPENDIXA - -- - _ _- - - -- - 135
APPENDIX B _- - -- - - - 146
REFERENCES. _ ..... -- -- _ __-..-_- - 162

LIST OF TABLES

Table 2.]: Values for Eq. (2.23) ....................................................................................... 14
Table 2.2: Face to node assignment for linear .................................................................. 15
Table 2.3: Values for Eq. (2.47) ....................................................................................... 21
Table 3.1: Engine geometry and operating conditions ..................................................... 55
Table 4.1: Parameterized piston model feature labels ...................................................... 60
Table 4.2: Heat transfer coefficient and ambient temperature values for parameterized
piston ............................................................................................................... 64
Table 4.3: Maximum and minimum temperature for the three sets of heat transfer
coefficients ...................................................................................................... 68
Table 4.4: Maximum and mean deformations on skirt ..................................................... 80
Table 4.5: Mesh and computation time comparison ......................................................... 80
Table 5.1: Flags for thermal analysis ................................................................................ 91
Table 5.2: Flags for static analysis .................................................................................... 92
Table 6.1: Engine geometry and operating conditions ................................................... 105
Table 6.2: Piston properties ............................................................................................ 106
Table 6.3: Material properties for aluminum alloy I345 ................................................ 107
Table 6.4: Thermal boundary conditions ........................................................................ 108
Table 6.5: Mesh sizes ...................................................................................................... 110

xi

LIST OF FIGURES

Figure 2.1: Piston CAD model, (a) top side, (b) underside ................................................ 6
Figure 2.2: Finite element mesh of piston, (a) top side, (b) underside ............................... 7
Figure 2.3: Linear triangular element ............................................................................... 1 1
Figure 2.4: Triangular element in its natural coordinates ................................................. 13
Figure 2.5: Linear tetrahedral element .............................................................................. 14
Figure 2.6: Nodal degrees of freedom for displacement analysis ..................................... 18
Figure 2.7: Tetrahedron natural coordinates ..................................................................... 19
Figure 2.8: 3D arbitrary body ........................................................................................... 22
Figure 2.9: Body and coordinate system for Reynolds lubrication equation .................... 38
Figure 2.10: 2D arbitrary body ......................................................................................... 39
Figure 2.11: Graphical representation of the population of a sparse matrix ..................... 46
Figure 2.12: Conjugate gradient method flowchart .......................................................... 49
Figure 3.1: Piston terminology ......................................................................................... 50
Figure 3.2: Piston assembly .............................................................................................. 52
Figure 3.3: Piston position extremities: bottom dead center and top dead center ............ 53
Figure 3.4: Piston position ................................................................................................ 55
Figure 3.5: Piston velocity ................................................................................................ 56
Figure 3.6: Piston acceleration .......................................................................................... 56
Figure 3.7: Angle (a ........................................................................................................... 57
Figure 4.1: Full piston model mesh in COSMOSDesignSTAR, (a) top side,

(b) underside ................................................................................................... 59
Figure 4.2: Parameterized piston model mesh in CASE, (a) top side, (b) underside ....... 59

xii

Figure 4.3: Parameterized piston dimensions ................................................................... 60
Figure 4.4: Pressure trace of optical engine ...................................................................... 61
Figure 4.5: In-cylinder air temperature for optical engine ................................................ 62

Figure 4.6: Heat transfer coefficient above piston crown using Annand and Woschni
correlations ...................................................................................................... 62

Figure 4.7: Heat transfer coefficients for parameterized piston ....................................... 64

Figure 4.8: Full piston model temperature distribution with boundary conditions of
Table 4.2 (temperature in Kelvin) ................................................................... 65

Figure 4.9: Full piston model temperature distribution with heat transfer coefficient
values increased by 15% (temperature in Kelvin) .......................................... 65

Figure 4.10: Full piston model temperature distribution with heat transfer coefficient
values decreased by 15% (temperature in Kelvin) ......................................... 66

Figure 4.1 l: Parameterized piston model temperature distribution with boundary
conditions of Table 4.2 (temperature in Kelvin) ............................................. 66

Figure 4.12: Parameterized piston model temperature distribution with heat transfer
coefficient values increased by 15% (temperature in Kelvin) ........................ 67

Figure 4.13: Parameterized piston model temperature distribution with heat transfer
coefficient values decreased by 15% (temperature in Kelvin) ....................... 67

Figure 4.14: Full piston model deformation due to thermal loading
(URES: resultant deformation in mm, deformation scale 1:50) ..................... 69

Figure 4.15: Parameterized piston model deformation due to thermal loading
(URES: resultant deformation in mm, deformation scale 1:50) ..................... 69

Figure 4.16: Mechanical loading, red arrows: atmospheric pressure, green arrows:
constraint, body load not shown ..................................................................... 70

Figure 4.17: Full piston model deformation due to mechanical loading
(URES: resultant deformation in mm, deformation scale 1:50) ..................... 71

Figure 4.18: Parameterized piston model deformation due to mechanical loading

(URES: resultant deformation in mm, deformation scale 1:50) ..................... 71
Figure 4.19: Pressure loads on skirt .................................................................................. 73
Figure 4.20: Full piston model deformation due to load case 1 ........................................ 73
Figure 4.21: Parameterized piston model deformation due to load case 1 ....................... 74

xiii

Figure 4.22: Full piston model deformation due to load case 2 ........................................ 74

Figure 4.23: Parameterized piston model deformation due to load case 2 ....................... 75
Figure 4.24: Full piston model deformation due to load case 3 ........................................ 75
Figure 4.25: Parameterized piston model deformation due to load case 3 ....................... 76
Figure 4.26: Full piston model deformation due to load case 4 ........................................ 76
Figure 4.27: Parameterized piston model deformation due to load case 4 ....................... 77
Figure 4.28: Full piston model deformation due to load case 5 ........................................ 77
Figure 4.29: Parameterized piston model deformation due to load case 5 ....................... 78
Figure 4.30: Full piston model deformation due to load case 6 ........................................ 78
Figure 4.31: Parameterized piston model deformation due to load case 6 ....................... 79

Figure 4.32: Parameterized piston model temperature distribution from CASE
(temperature in degrees Celsius) ..................................................................... 81

Figure 4.33: Parameterized piston model deformation obtained from CASE

(deformation in mm) ....................................................................................... 82
Figure 4.34: Skirt oil film thickness on thrust side at 33 CAD ........................................ 83
Figure 4.35: Skirt oil film thickness on anti-thrust side at 33 CAD ................................. 84
Figure 4.36: Friction force between piston and cylinder liner .......................................... 84
Figure 4.37: Contact points of piston with cylinder wall .................................................. 85
Figure 4.38: Cyclic piston rotation about pin ................................................................... 86
Figure 4.39: Cyclic piston pin lateral motion ................................................................... 86
Figure 4.40: Experimental oil film thickness measurement at TDC of intake stroke ....... 88
Figure 4.41: Simulated oil film thickness prediction at TDC of intake stroke ................. 88
Figure 5.1: Piston orientation ............................................................................................ 90
Figure 5.2: Overview of the 1t-fea program flow .............................................................. 93
Figure 5.3: 2D skirt profile ............................................................................................... 95
Figure 5.4: Forces and moments acting on piston ............................................................ 98

xiv

Figure 5.5: Nodal temperature as predicted by n-fea and COSMOSDesignSTAR ........ 104

Figure 5.6: Nodal resultant deformation as predicted by rt-fea and COSMOSDesignSTAR

....................................................................................................................... 104
Figure 6.1: Pressure trace used for simulation ................................................................ 106
Figure 6.2: Piston CAD model, (a) isometric view, anti-thrust side shown,

(b) underside, (c) front view, anti-thrust side shown, ((1) side View,

red: thrust side ............................................................................................... 107
Figure 6.3: Imported piston meshed geometry, global element size: 9.2502 mm, skirt

element size: 3 mm ....................................................................................... 109
Figure 6.4: Skirt mesh, (a) thrust side, (b) anti-thrust side ............................................. 109
Figure 6.5: Hydrodynamic force on thrust side for meshes of Table 6.5 ....................... 1 10
Figure 6.6: Hydrodynamic force on anti-thrust side for meshes of Table 6.5 ................ 1 11
Figure 6.7: Total hydrodynamic fiiction force for meshes of Table 6.5 ......................... 111
Figure 6.8: Hydrodynamic force over three cycles, C = 100 um ................................... 112
Figure 6.9: Hydrodynamic force over three cycles, C = 50 um ..................................... 113
Figure 6.10: Hydrodynamic pressure convergence at 31rd cycle, C = 100 um ................. 113
Figure 6.11: Hydrodynamic pressure convergence at 3rd cycle, C = 50 um ................... 114
Figure 6.12: Temperature distribution ............................................................................ 115

Figure 6.13: Deformation at 405 deg.: x-component, (a) C = 100 mm, (b) C = 50 um .. 116

Figure 6.14: Equivalent strain at 405 deg., C = 100 um ................................................. 116
Figure 6.15: von Mises stress at 405 deg., C = 100 um .................................................. 117
Figure 6.16: 3D von Mises yield criterion at 405 deg., C = 100 um .............................. 118
Figure 6.17: Hydrodynamic pressure on thrust side at 315 deg., C = 100 um ............... 1 19

Figure 6.18: Hydrodynamic pressure on anti-thrust side at 315 deg., C = 100 um ........ 1 19
Figure 6.19: Hydrodynamic pressure on thrust side at 405 deg., C = 100 um ............... 120
Figure 6.20: Hydrodynamic pressure on anti-thrust side at 405 deg., C = 100 um ........ 120

Figure 6.21: Hydrodynamic pressure on thrust side at 405 deg., C = 50 um ................. 121

XV

Figure 6.22: Hydrodynamic pressure on thrust side at 405 deg, C = 50 pm ................. 121

Figure 6.23: Oil film thickness on thrust side at 405 deg., C = 100 um ......................... 122
Figure 6.24: Oil film thickness on anti-thrust side at 405 deg., C = 100 um .................. 123
Figure 6.25: Total hydrodynamic force on piston .......................................................... 124
Figure 6.26: Plane parallel to pinhole axis and coincident with piston axis ................... 124

Figure 6.27: Histogram of resultant deformation on thrust and anti-thrust side due to

thermal expansion only ................................................................................. 125
Figure 6.28: Total contact force on piston ...................................................................... 126
Figure 6.29: Total friction force on piston ...................................................................... 126
Figure 6.30: Total moment on piston .............................................................................. 127
Figure A.1: Flow chart for the main program of rt-fea (Part 1 of 2) ............................... 135
Figure A.2: Flow chart for the main program of rt-fea (Part 2 of 2) ............................... 136
Figure A.3: Flowchart for FEMmeshread subroutine ..................................................... 137
Figure A.4: Flowchart for FEMthermal subroutine ........................................................ 138
Figure A.5: Flowchart for FEMstatic subroutine (Part 1 of 2) ....................................... 139
Figure A.6: Flowchart for FEMstatic subroutine (Part 2 of 2) ....................................... 140
Figure A.7: Flowchart for FEMGuyan subroutine ......................................................... 141
Figure A.8: Flowchart for FEMstaticL UB subroutine .................................................... 142
Figure A.9: Flowchart for ReEq subroutine .................................................................... 143
Figure A.10: Flowchart for FEMstrain subroutine ......................................................... 144
Figure B]: Temperature distribution ............................................................................. 146
Figure B.2: Deformation: x-component .......................................................................... 147
Figure B.3: Deformation: y-component .......................................................................... 147
Figure B.4: Deformation: z-component .......................................................................... 148
Figure B.5: Deformation: resultant ................................................................................. 148
Figure 8.6: Deformation: resultant, scale 1:800 ............................................................. 149

xvi

Figure B.7: Normal strain: x-component ........................................................................ 149

Figure B.8: Normal strain: y-component ........................................................................ 150
Figure 89: Normal strain: z-component ........................................................................ 150
Figure 8.10: Shear strain: xy-component ........................................................................ 151
Figure 8.11: Shear strain: xz-component ........................................................................ 151
Figure 3.12: Shear strain: yz—component ........................................................................ 152
Figure 8.13: Equivalent strain ........................................................................................ 152
Figure 8.14: Normal stress: x-component ...................................................................... 153
Figure B.15: Normal stress: y-component ...................................................................... 153
Figure B.16: Normal stress: z-component ...................................................................... 154
Figure B.17: Shear stress: xy-component ........................................................................ 154
Figure B. 1 8: Shear stress: xz-component ........................................................................ 155
Figure B.19: Shear stress: yz-component ........................................................................ 155
Figure 8.20: Principle stress: 1St principal direction ....................................................... 156
Figure B.21: Principle stress: 2'"d principal direction ...................................................... 156
Figure B.22: Principle stress: 3rd principal direction ...................................................... 157
Figure B.23: von Mises stress ......................................................................................... 157
Figure 8.24: Stress intensity ........................................................................................... 158
Figure B.25: 3D von Mises yield criterion ..................................................................... 158
Figure B.26: Hydrodynamic pressure: thrust side .......................................................... 159
Figure 8.27: Hydrodynamic pressure: anti-thrust side ................................................... 159
Figure 8.28: Oil film thickness: thrust side .................................................................... 160
Figure 8.29: Oil film thickness: anti-thrust side ............................................................. 160

xvii

2D

3D

A

(11, 02. a3

01. a2, a3, a4

BDC

b:

C1» 02. C3

CI, Cz. C3. C4
C

CAD

CAE
CASE

const

KEY TO SYMBOLS AND ABBREVIATIONS

gradient, divergence

two dimensional

three dimensional

element area

shape function constant for triangular element
shape function constant for tetrahedral element
piston acceleration

skirt area

shape function x-coefficient for triangular element

shape function x-coefficient for tetrahedral element

matrix of shape function derivatives

bottom dead center

body load

shape function y-coefficient for triangular element

shape function y-coefficient for tetrahedral element

nominal piston to cylinder bore clearance, compliance matrix
Computer Aided Design

Computer Aided Engineering

Cylinder-kit Analysis System for Engines

constant

pin offset, total deformation

xviii

db

dim

FEA

FEM

F/‘c

F/h

piston deformation due to body load

piston deformation due to unit body load

degrees

piston deformation due to combustion gas pressure
piston deformation due to unit pressure on crown
resultant deformation

total deformation on skirt side m

deformation on skirt due to hydrodynamic pressure loading
deformation due to thermal strain

modulus of elasticity

function

function

contact force

Finite Element Analysis

Finite Element Method

total fiiction force

contact fiiction force

hydrodynamic friction force

hydrodynamic force

total force

load vector due to physical wedge

load vector of element ‘e’ due to squeeze film

xix

kx, ky, kz

[Kred]

L1. L2. L3
L1. L2. L3. L4
m

M

load vector due to squeeze film

load vector of element e due to physical wedge

heat transfer coefficient, oil film thickness
heat transfer coefficient of piston ring-pack
heat transfer coefficient of piston skirt

heat transfer coefficient of piston underside

combustion gas heat transfer coefficient

squeeze film velocity
functional
1'“ node

1"” node

global stiffness matrix
stiffness matrix of element e

thermal conductivity vector
thermal conductivity

thermal conductivities along x, y, z axes

reduced global stiffness matrix

connecting rod length

natural coordinates for triangular element
natural coordinates for tetrahedral element
number of matrix entries

total moment

XX

min

N

[N]

n

N1. N2. N3

N1. N2. N3. N4
nel

nx, n,, it;

P

p<e>

P

Pres

contact moment

moment due to contact friction force
moment due to hydrodynamic friction
hydrodynamic moment

minute

engine speed

shape function matrix

index number

shape functions for triangular element
shape functions for tetrahedral element
number of elements

direction cosines

pressure

load vector of element 6

load vector, pressure vector

load vector of element 6 due to body forces

load vector due to body forces

contact pressure
combustion gas pressure
hydrodynamic pressure

pressure residual

load vector of element e due to surface forces

xxi

(Qt

[R]

rev

subscript m
superscript e
To.

TI

T 2

T 3

T

T09)

avg

ref

TDC

load vector due to surface forces

load vector of element 6 due to thermal strains

load vector due to thermal strains

total pressure

uniform pressure

vector of nodal displacements
elasticity matrix

residual

residual, crankshaft radius
revolution

search direction

skirt side m

element e

ambient temperature

ambient temperature of piston ring-pack
ambient temperature of piston skirt
ambient temperature of piston underside

temperature

average temperature of element ‘e’
reference temperature

top dead center

combustion gas temperature

xxii

displacement vector

components of displacement parallel to x, y, z axes
element volume

relative velocity

velocity of body a

velocity of body b

piston velocity

x-axis, x-direction, x-coordinate

local skirt x—coordinate

variable, y-axis, y-direction, y-coordinate
piston position

local skirt y-coordinate

z-axis, z-direction, z-coordinate

xxiii

an r

a, 3, y, 6

PO

01. 02. 03

Gin!

coefficient of thermal expansion
shape function power for triangular element
shape function power for tetrahedral element
constant for conjugate gradient
shear strain

surface boundary

wavy deformation

change in temperature

time step

normal strain

equivalent strain

initial thermal strain

1'” generalized coefficient
engine crank angle

step length

strain energy

lubricant viscosity

Poisson’s ratio

variable

potential energy

normal stress

principal stresses

stress intensity

xxiv

01'0"

“Na

'9' 'S-II

“S

(Ur

von Mises stress
shear stress

angle of connecting rod with piston axis, body force

prescribed boundary condition

surface force vector

element behavior

switch function

work done on a body by external forces
weight function, angular speed
domain, waviness

underrelaxation factor

XXV

CHAPTER 1

INTRODUCTION

1.1 Motivation

As the oil prices increase, the emission standards tighten, and the competition in
the automotive market strengthens, the search for a more efficient reciprocating internal
combustion engine becomes vital. The prime “navigators” in this search are
computational tools. Such tools allow for a fast and relatively cheap course for a
prototype design, but can also be utilized for troubleshooting and optimization of existing
designs.

The internal combustion engine has been around for more than a century but its
physics are still not well understood. Its cyclic behavior represents an unsteady complex
system with multiple physical processes occurring simultaneously. The system then is a
multidisciplinary one which can be described by combustion [21, 41], thermodynamics
[3, 39], heat transfer [27, 29], solid mechanics [2, 53], fluid mechanics [43, 46], dynamics
[22, 26], and tribology [32, 38].

The “heart” of a reciprocating internal combustion engine is the piston. The
function of the piston is to convert the thermal energy of the combustion gases into the
mechanical energy that drives the engine. It is generally believed that about half of the
mechanical energy is lost just at the piston assembly. Therefore piston design is a very
important factor in engine efficiency.

The piston is confined within a cylinder bore, which acts as a guide to the piston’s

reciprocating motion. The piston skirt slides on the cylinder’s surface, thus providing the

support for the whole piston body. Under favorable operating conditions an oil film
separates the two surfaces. During a cycle the piston is exposed to sharp temperature and
pressure gradients due to the combustion gases just above its crown. The reciprocating
motion also exposes the piston to inertial loads. These gradients and loads affect the
mechanical behavior of the piston skirt, which directly affects the thickness of the oil
film, which in turn, determines the piston dynamics and the interaction of the piston with
the cylinder wall. If the oil film thickness is too low, then piston scuffing can be a
problem. If the thickness is too high, piston slap can be a problem. Therefore piston
quality depends on optimal geometrical design, body mass and material selection, that
ensure minimal energy loss due to inertia, heat transfer and piston-bore interaction.
However finding the optimal set of design characteristics that would yield optimal piston
quality is a very challenging job. Consequently computational tools become very

important in piston design.

1.2 Previous Efforts in Piston Modeling

Piston modeling can be divided into four subcategories: thermal modeling, solid
modeling, elastohydrodynamic lubrication modeling, and dynamic modeling. Over the
years several research efforts have been made in these areas.

Armand [1] and Woschni [55] proposed correlations to calculate the heat transfer
coefficient of the combustion gases above the piston crown. Woschni [56] did some
further work for the evaluation of the heat transfer coefficients for a high speed diesel
engine piston. He used the electrolytic tank analogue to obtain temperature measurements

to calculate the heat transfer coefficients. He proposed ranges for heat transfer

coefficients describing the different surfaces of the piston. Wu et al. [57] modified
Annand’s correlation to include radiation effects, and they developed a numerical model
for the calculation of the temperature distribution. Their analytical results showed
reasonable agreement with experimental ones.

Li [35] considered the piston’s thermoelastic behavior. He assumed that
temperature fluctuations during a cycle affect a piston layer only about 2 mm thick.
Beyond this layer the temperature is steady, given enough operating time for the engine.
Consequently he treated the piston’s thermoelastic behavior as a steady state problem. He
used experimental temperature measurements to propose a range of heat transfer
coefficients for an aluminum gasoline piston. He used these results and a finite element
model of a quarter of the piston to investigate the thermoelastic behavior of the piston
over a cycle.

Li et a1. [37] developed an automotive piston lubrication model to study the
effects of piston pin location, piston-to-cylinder clearances and lubricant viscosities on
piston dynamics and friction. They solved for particular solutions of the Reynolds
equation using finite differences, and used the Newton-Raphson method to solve for the
nonlinear equations of motion for the piston. In their model they assumed a rigid piston
model. Li [36] considered the elastic deformation of the piston skirt; integrating this with
hydrodynamic lubrication has formed the elastohydrodynamic lubrication analysis which
is considered by most of the recent efforts [9, 10, 11, 18, 30, 31, 40, 45, 54, 58, 59]. All
these have contributed in solving the elastic deformation of the piston using the finite
element method. Oh et al. [40] used the finite element method to solve for the Reynolds

equation; however they linearized it, thus solving for a set of linear equations. The rest of

 

the efforts use finite differences to solve for the Reynolds equation. This method however
requires mapping of nodal information, back and forth, from the finite element mesh to
the finite difference grid. This can result in the loss of crucial numerical information,
especially where sharp gradients exist between two nodes. Zhu et al. [58, 59] were the
first to consider the elastic deformation of the cylinder bore. They also proposed a
formula for the asperity contact pressure of an aluminum piston against a cast iron
cylinder according to Johnson’s model (1985). They also used the Reynolds equation
developed by Patir and Cheng (1978) which accounts for the effect of surface texture on
hydrodynamic lubrication. Duyar et al. [11] coupled the Reynolds equation with the
mass-conserving Reynolds equation to solve for the hydrodynamic pressure using finite
volumes. Computer aided engineering (CAE) tools for piston analysis have been

introduced in recent years, such as CASE [4], PISDYN [l 1] and GLIDE [24].

1.3 Objective

In this thesis, the author joins the disciplines of heat transfer, solid mechanics,
fluid mechanics and tribology to develop a computational model for the assessment of
piston quality. The Cylinder-kit Analysis System for Engines (CASE) is a comprehensive
cylinder-kit simulation software that is being widely used to predict piston dynamics,
ring-pack dynamics, oil film thickness and friction [6, 7, 12, l3, 14, 15, 16, 42]. Currently
CASE uses a parameterized piston model which is generated using a set of selected
geometrical parameters from the full piston model. The main advantage of the
parameterized model is its relatively simple geometry. This allows for the construction of

a simple small mesh where the governing equations can be solved very fast. Another

 

advantage of the parameterized model is its self-generation. This allows CASE, when
coupled with optimization sofiware, to perform an optimization simulation very quickly
without the need of external computer-aided design (CAD) software to change the
piston’s geometry. The author develops a piston finite element analysis program to
supplement the existing one of CASE. This program is capable of importing the meshed
geometry of a full CAD piston model and performs a cyclic analysis to assess the piston
quality. The parameterized approach has proven to be a less expensive approach in piston
quality assessment in terms of cost, time and effort. The objective then is to be able to
assess the optimal piston design characteristics like piston-pin location, skirt profile and
piston ovality from the parameterized piston model and use them to update the CAD full
piston model. Then a simulation for the modified CAD full piston model can be
conducted to verify if the selection of optimal design characteristics can improve the
performance.

In this thesis only the finite element part of the program is presented. At this stage
piston secondary dynamics are ignored. The piston is assumed to be moving right in the
centre of the cylinder bore, with no transverse movement or tilting. Viscosity variations

along the length of the cylinder bore and cylinder bore deformation are ignored.

 

CHAPTER 2

THE FINITE ELEMENT METHOD

2.1 Introduction

This chapter introduces the reader to the finite element method, from hereon
referred to as FEM, used in the preparation of this work. The philosophy of FEM is to
find the solution of a complicated problem by replacing it by a simpler one. Consider the
piston shown in Figure 2.1; it has a very complex geometry. The existing mathematical
tools are not sufficient to find an exact solution for its behavior (temperature distribution,

displacement distribution) under thermal and mechanical loading.

 

(a) (b)

Figure 2.1: Piston CAD model, (a) top side, (b) underside

However an approximate solution can be obtained using FEM. The complex
geometry of the piston can be divided into an assembly of many small sub-geometries

called finite elements (Figure 2.2). These elements are interconnected at specified joints

called nodes or nodal points. These nodes are usually found on the element boundaries
where adjacent elements are connected. Again the behavior of these finite elements is not
known, however it is assumed that it can be approximated by simple functions. These
approximating functions, also known as shape functions, are defined by the element
behavior at the nodes. Equilibrium equations can be readily developed to describe the
behavior of these finite elements. The unknowns arising from the development of these
equations are the nodal values of the element behavior. The simultaneous solution of
these equations describes the approximate behavior of the whole geometry. One
advantage of FEM is that the approximate solution can be improved by refining the mesh.
Of course however this would require more computational time. All of the theories
introduced in this chapter are extracted from literature [5, 20, 23, 28, 32, 33, 34, 38, 44,

47, 48, 49, 50, 51, 52,].

 

(a) (b)

Figure 2.2: Finite element mesh of piston, (a) top side, (b) underside

2.2 Approximate Solutions
Consider the function given in Eq. (2.1)
f(0) = y"(Q)+y'(Q)+y(Q) (2.1)
The residual, R, is given by
R = y"(Q) +y'(Q)+y(Q) -f(Q) (2.2)
The exact solution of Eq. (2.1) would yield R = 0. With approximate solutions, if

R cannot be made zero, it is brought as close to zero as possible.

2.2.1 Galerkin’s Method

This method, also known as the method of weighted residuals, tries to reduce the
residual over the entire domain by making its weighted average zero with respect to as
many independent weighting functions, a), as there are unknown parameters. Assume that

Eq. (2.1) is approximated by n unknown parameters, 3. Thus,
R=R(Q, 50, 51, ..., E") (2.3)
This results in,

MR“), 5,, 52, 5,) d.(2=0

.0
[@Rm, 5,, 52, 5,) d.(2=0

{2 (2.4)

[(0, R(.(2, 5,, 52, 5,) d0=0
.(2
Galerkin’s method uses the element shape functions as the weighting functions.

This ensures that the weighting functions are independent of each other and consequently

the resulting algebraic equations are independent.

2.2.2 Variational Method

This method arises from calculus of variations. The main idea behind calculus of
variations is to determine the extreme (maximum and minimum) or stationary values of
functionals. A functional, 1, is defined as function of other functions. Functionals of
equilibrium equations have been derived and are found in literature. Considering Eq.

(2.1) again approximated by n unknown parameters, 5', its functional can be given by,

[2 ]F(9, 50, 5,, E,)d.(2 (2.5)
.0

Values of 5,, for i = l, 2,. . ., n, are sought which will minimize I, that is,

 

ar ar 31
—=—=...= =0 2.6
as, 35, a5, ( )

Equation (2.6) will result in n equations with n unknowns.

2.3 Elements and their Shape Functions
2.3.1 Linear Triangular Element

A linear triangular element (Figure 2.3) is a three-node element with straight
edges. The nodes are numbered in a counter-clockwise order. Node 1 can be specified
arbitrarily. Let the element behavior (temperature, displacement) be CD. Inside the element
(D is assumed to vary linearly, Eq. (2.7), and 45,, (Dz, $3, Eq. (2.8), denote the behavior

value at the respective nodes.

¢(x,y)=rh +222 mm (2.7)

Q=m+ma+mn
$2 = 7714-772 1'2 +773y2 (2.8)
4’3 = 771+772 x3 “13%

The solution of Eqs. (2.8) leads to
1
771 = 31(014’1 +(112452 + “3%)

l
772 = 350291 +b2¢2 +b3¢3) (2.9)

1
'73 = 5(614’1 +62% ”3453)

where A is the area of triangle 1-2-3. The area is given by

l
1 x2 y2 (2-10)
1

and

m=en-en

%=&M‘nh

%=Mh-hh

fi=n-h

la=yy<n (21w
%=n—h

q=e-n

Q=M-B

C3 =X2 —x,

m

 

 

Figure 2.3: Linear triangular element

Substituting Eq. (2.9) into Eq. (2.7) and rearranging leads to
¢(x,y)= N1(x,y)¢1 + N2(x.y)¢2 + N3(x.y)¢3

where,

N] =312(a,+b1x+c,y)

N2 =2—1/I(a2 +b2x+c2y)

N3 =3!;1'(03 +b3x+c3y)

(2.12)

(2.13)

N ,, N 2, and N3 are the element shape functions. It can be shown that evaluating the

shape functions at x,, y], N 1 = 1, whereas N1 = 0, at nodes 2 and 3. Similarly the shape

function N2 and N 3 have a value of 1 at nodes 2 and 3 respectively and 0 at the other

nodes.

11

Now let,

[N]=[N, N2 N3] (2.14)
and
qb,
45(9): Q (2.15)
4’3

where, OM is the vector of nodal unknowns of element e. Thus,

ab = [N]O(e) (2.16)

Since the functions are linear, their derivative with respect to x or y is constant.
This suggests a constant derivative of O within the element, Eq. (2.17). Consequently
very small elements have to be used where sharp gradients are expected in the value of O.

aab —?__N1_¢, a__1v2 a, +_a__~,

O 2.17
T): 8x 3x 3 ( )

Now consider a triangular element in its natural coordinates (Figure 2.4), L 1. L 2,

L 3. These three coordinates define a point P in terms of the areas, Eq. (2.18),

A A A
=.Ai,1,2=_AZ.,1,3=_/.13. (2.18)
and
A,+A2+A3=A (2.19)
thus,
L,+l/2+L3=l (2.20)

12

 

Figure 2.4: Triangular element in its natural coordinates

It can readily be shown that
L1=N],1/2=N2,L3=N3 (2.21)
The Cartesian coordinates of P are related to the natural ones by,

x = x1L1+ x2112 + x313

 

 

(2.22)
y = y1L1+Y2Q + y3l/3
For integrating polynomial terms in area coordinates Eq. (2.23) can be used,
I I 1
[[Lf’nggdAz “#7" 2A (2.23)
A (61+ ,6 + y+ 2)!
and since Eq. (2.21) is true, then
I 1 I
[[NFNfod/r: “'6'” 2A (2.24)
(CH ,6 + 7+ 2)!

A

13

Table 2.1 gives some common values for the integral of Eq. (2.23).

Table 2.1: Values for Eq. (2.23)

 

 

Value of Value of integral in Eq. (2.23)/A

a A 7

0 0 0 1

l 0 0 1/3
2 0 0 1/6
1 1 0 1/12
3 O 0 l/ 10
2 1 0 1/30
1 l 1 1/60
4 0 0 1/15
3 , 0 1/60

 

 

 

 

2.3.2 Linear Tetrahedral Element

The linear tetrahedron element is made of four flat triangular elements or faces,
interconnected at the four comers, nodes 1, 2, 3 and 4. The nodes on any face are
numbered in a counter—clockwise manner as seen from the vertex opposite the face. Table

2.2 shows the faces with their respective nodes as used in this work.

YA

 

—.L

Figure 2.5: Linear tetrahedral element

14

Table 2.2: Face to node assignment for linear

tetrahedral element

 

 

Face Nodes
1 1-2-3
2 1-4-2
3 2-4-3
4 1-3-4

 

 

The derivation of the element shape firnctions is very similar as for the triangular
element. Again here the element behavior, O, is assumed to vary linearly within the

element,
¢(x,y,2)=771+772x+773y+774z (2.25)

@=m+ma+mn+ma
4’2 = 771 +772 x2 +773)? +0422

(2.26)
4’3 = ’71 +772 1‘3 +773J’3 +0423
4’4 = 771 + 772 x4 + 773M + 77424
The solution of Eqs. (2.26) leads to
l
'71 = 37(014’1 +42% +a34’3 + 044’4)
l
772 = —(b1¢1 +bz¢2 +b39153 +b4¢4l
2V (2.27)

l
773 = 2—V(Cl¢1 +Cz¢2 +C3¢3 +C4¢4)

l
774 =§7ld1¢1 +“’2952 +43% “14454)

15

where V is the volume of the tetrahedron 1-2-3-4. The volume is given by,

1‘1
x2

x3

p—a
y—ap—‘u—‘p—s

 

x4

and,

Y1
Y2
Y3
Y4

Y2
Y3
Y4

21
22
Z3
Z4

 

22

24

1
l
l

(2.28)

(2.29)

(2.30)

(2.31)

(2.32)

for i = 1,2, ..., 4. The subscripts of x, y, and z, in Eqs. (2.29) to (2.32), are interchanged

cyclically in the order 4, l, 2, 3. These equations are valid for a tetrahedral numbered as

in Figure 2.5.

Substituting Eqs. (2.27) into Eq. (2.25) leads to,

¢(X,Y,Z)= Nl(X,Y9Z)¢1 +N2(anaZ)¢2 +N3(st9z)¢3 +N4(staZ)¢4 (233)

where,

16

N1=—l—(a1+b, x+cl y+d1 2)

6V
N2 =-1—'(02 +b2x+C2y+d22)

61V (2.34)
N3 =—(a3+b3x+c3y+d3z)

6V
N4 =817(a4+b4x+c4y+d4z)

N I, N 2, N 3, and N4 have a value of one at nodes 1, 2, 3, and 4 respectively and zero

 

 

elsewhere.
Now letting
[N]=[N, N2 N3 N4] (2.35)
and
I¢1‘
_ ab
db“): 2) (2.36)
4’3
IP41
then
(D = [N]é(") (2.37)

As in the case of the triangular element the derivative of O with respect to x, y or
2, within the element is constant, Eq. (2.38).

a_a> _8_N,¢1+ 6__N2¢2 +8_N3¢3+E_9_N_4

(D 2.38
8x 8x 8x 8x 4 ( )

So far O is assumed to be a scalar quantity; however there are cases (displacement

analysis) where O is a vector (Figure 2.6).

17

 

Figure 2.6: Nodal degrees of freedom for displacement analysis

In such a case O(e) becomes,

65(3) =< v2 ) (2.39)

 

 

and,

18

 

 

[N1 0 0
N2 0 0
N3 0 0
N4 0 0
0 N, 0

[N]T = O N’- O (2.40)
0 N3 0
0 N4 0
0 0 N,
0 0 N2
0 0 N3
[0 0 N4_

so that,

u

45: v =[N]O(e) (2.41)
W

Like the triangular element a tetrahedron can be described in its natural or volume

coordinates, L 1, L2, L 3, L4. These coordinates define a point P in terms of volumes (Figure

2.7).

 

Figure 2.7: Tetrahedron natural coordinates

l9

where,
V= Volume of 1-2-3-4
V, = Volume of P-2-3-4
V2 = Volume of P-1-3-4
V3 = Volume of P-1-2-4
V4 = Volume of P-1-2-3
thus,
fi+fi+fi+fl=V
and
Q+Q+Q+Q=l
It can be readily shown that
L1 =1V1 2 12=N2a 113=N32 L4=N4
The Cartesian coordinates of P are related to the natural ones by,

x = x1L1+x21Q +X3L3 +X4L4
y = y1L1 + y21/2 + y31/3 + Y4L4
Z 1' 21L] '1' Zzbz + 231/3 '1' 241.4

(2.42)

(2.43)

(2.44)

(2.45)

(2.46)

For integrating polynomial terms in volume coordinates Eq.(2.47) can be used,

 

 

. . .6!
La 3 7146de a'fl'y'
911L213 4 (a+fl+y+6+3)!
and since Eq. (2.45) is true, then
1 1 151
NaNflN7N5 dV= a'fl‘y' '
{III 2 3 4 (a+fl+y+6+3)!

20

(2.47)

(2.48)

Table 2.3 gives some common values for the integral of Eq. (2.47)

Table 2.3: Values for Eq. (2.47)

Value of Value of integral in Eq.
(2.47)/A
1
1/4

1/10

1/20

1/20

1/60

1/ 120

1/35

 

 

bt—Nw—Nr—Oh
o-——‘o~—c>ccata
o-—<:c>c>o<:c:w
o<:c>c>o<:c>c>m

 

 

 

 

 

2.4 Finite Element Formulation
2.4.1 Heat Equation
For a three-dimensional steady state problem with no heat generation (Figure 2.8)

the governing differential equation is,

V(IEVT) =0 (2.49)

with a convective boundary,

3T 3T 3T
kxgx—nx+ky-é—y—ny+kz-é:nz+h(T—T°o)=0 OIIF (2.50)

21

 

Figure 2.8: 3D arbitrary body

For an isotropic material where k, = k, = k2 = k, (2.49) reduces to,

V2T=O
and (2.50) to,
8T 3T +8T
an, +3y 11+), 82 —nz +— (—T—T°°)=0

(2.51)

(2.52)

Using Galerkin’s method, multiplying by the weight function and integrating over the

domain, Eq. (2.51) becomes

[ szT do“) = 0
0(9)

Using Green’s first identity,

[(VT Vw+a1V2T T=)d.(2 ][—n +%§ny+g—Tn )wdl‘
.0

Eq. (2.53) yields,

] VTdeQ‘e) = [

[LT at" 87"
0(9) [(6’)

(9)
8x "x x+a—ny+a—nz ,jwdl‘

22

(2.53)

(2.54)

(2.55)

From Eq. (2.52), Eq. (2.55) becomes,

[ VT-Va) 2112(8) =— [ fl(T—T,,) wdl‘(e) (2.56)
{2(9) r(9) k
Rearranging yields,
[ VT-Va) 210+}: [ Ta) dl‘zfl- ] wdrle) (2.57)
(e) k (e) k (e)
.{2 1" F
By approximating T,
T 2 N,- 7} (2.58)
and,
a: . 2,. in
8x Bx
a: = I} 8_N,_ (2.59)
8y 8y
8T 8N~
_ z T- __L
32 l 82
and,
a). = N. (2.60)

Eq. (2.57) through (2.60) yield,

 

h
d.Q(e) +— N-N- arr“) Tie)
3x 3x By By 82 82 k I I J {' }

”(e1 rte)

(2.61)

Equation (2.61) can be expressed as

([Kfe)]+[1<§e)]) Tie) = 13“?) (2.62)

23

with,

 

[KM]: [ aN" aNj+aNi aNj+aNi 8N1 dale) (2.63)
1 0(8) 3x 8x 8y By 82 dz
(1 _h
[5129}; [ N,Nj dF(") (2.64)
fl?)
13(8):.i-Tr- [ N, dl‘(e) (2.65)
k rlt’)

Now recalling the shape functions for a tetrahedral element, Eq. (2.34), the shape

function derivatives are given by,

.a_.NL -_—. i (2.66)
3x 6V
% = i (2.67)
By 6V
8N,- d,
__ = _ 2.68
dz 6V ( )
Thus considering the above equations, Eq. (2.63) becomes
- 1
[Kl(8):l= fl m(blbj +Cicj +d1d1) dg(e)
[2(8)
(2.69)
l
=6—[Y-(bibj +Cicj +didj)

Now considering Eq. (2.48) and assuming convection on face 1 of the tetrahedron

Eq. (2.64) becomes,

24

j
Flt)
2 1 1 0 (2.70)
h123A1231 2 1 0
l2k l l 2 O
_0 0 o 0‘

 

 

where h 2 2 3 is the heat transfer coefficient on face 1 and A 123 the area of face 1. The fourth
row and column are all zeros as M, is zero on face 1. There are three other forms of Eq.
(2.70) corresponding to faces 2, 3, and 4. Similarly the diagonal terms will be 2 and off-
diagonal terms 1. The terms of the row and column associated with the node not lying on
the face will be zero.

In a similar manner, for convection on face 1, Eq. (2.65) can be expressed as

13(9) :11: J. Ni dr(e)
k r(€)

(2.71)
= h123Too A123, 1,
3k 1

 

 

Again there are three other forms of Eq. (2.71) corresponding to each face. Here the row
corresponding to the node not lying on the face will be zero.

Eq. (2.69) to (2.71) represent the element matrices and vectors; assemblage of

these leads to the global stiffness matrix, [K], and global load vector P

”#:146114111

_ nel _.
P = Z P“) (2.73)

e=1

25

thus,

[K]T = 13 (2.74)

with T being the vector of the unknown nodal temperatures.

2.4.2 Static Analysis
In this section the general equations governing solid and structural mechanics are
presented. The principle of minimum potential energy, also known as displacement

method, is used in deriving the finite element equations.

Equilibrium Eflations

Consider the arbitrary body of Figure 2.8. If a load is applied on surface, F, and if
an element of material is considered within the domain [2, it must be in equilibrium due
to the internal stresses developed. This leads to the internal equilibrium equations for a

three-dimensional body,

a6” BQW+BQZ+

 

 

 

=0

dx + dy dz P
a: do a:

y" W W =0 2.75
dx + dy + dz +¢y ( )
asz 8121’ 80'22

=0
dx + dy + dz +¢z

where 0',“ and 0'22 are the normal stresses, 2'

x,,2',,1'_,z,2' 2’ and? are

, ayy , y] zx 9 yz a zy 9
the shear stresses, and ¢x , ¢y , and (9,, are the body forces per unit volume acting along

the x, y and 2 directions respectively. For an isotropic material,

26

1,0, = ryx
sz = sz (2.76)
Tyz = ’22

Thus for a three-dimensional isotropic problem there are only six independent terms,
0",“, 0'”, 022, 7,, , 7,2, and ryz.

Stress-strain Relations

Hooke’s law gives the stress-strain relations for a linearly elastic isotropic three-

dimensional solid as,
a=[R](a—a,) (2.77)

where 6 is the stress vector,

i (2.78)

Q1
11

 

 

where E is the strain vector,

7 " (2.79)

M
II

 

 

27

where 5‘0 is the initial strain stress vector given by a temperature differential AT , and the

coefficient of thermal expansion, a,

and the matrix [R] is given by,

where,

E is the Young’s modulus of elasticity, and v is the Poisson’s ratio.

 

 

 

 

 

 

 

 

 

8,0,0 [1
EYYo l
82 l
=< zo %=0'AT{ k
ny0 0
7x20 0
0

[7y20 L1
"A B B o 0 0“
B A B O 0 0
B B A O O 0
0 O O C O 0
O 0 O 0 C 0
_O 0 0 0 O C-
_ E<1-v)
—(1+V)(1—2v)

_ VE

(1+V)(1—2v)

E

C:

2(1+V)

28

(2.80)

(2.81)

(2.82)

Strain-displacement Relations

The deformation of any elastic body under any given loading can be described by

three displacement components, u, v, and w, parallel to the x, y, and 2 directions

respectively. The normal strains are given by,

du

gig—a?

dv

8y), = 5;

8-7

and the shear strains by,

_du
ny-ay

du
7x2

dv

=—+—
7y: dz dy

Boundngonditions

_dw

...Z 82

dv
dx
dw

=55;

dw

(2.83)

(2.84)

The boundary conditions can be either on displacements or on stresses. In this

work only boundary conditions on displacements are considered. Displacements of nodes

on boundary are set to zero.

Compatibility Equations

A body under a given loading should remain continuous, that is no cracks or gaps

should appear in the body. Consequently the deformations field should be continuous and

29

single-valued. This requires a definite relation between normal and shear strains, the

compatibility equations, Eq. (2.85).

 

82822 328.10: _ 327x:
2 2 ‘

dx dz dxdz (2.85)

1 a _[any_ 87... ay,,)=aze,,

 

 

 

 

 

 

2 dx dz dx dy ) dydz
ii 87’20) +37yz_ 37x2 _azgyv

2 dy dz dx dy ) dxdz
.1._ a _37xy+37yz +37... =£§£
2 dz dz dx dy dxdy

Principle of minimum potential energy

The principle of minimum potential energy states that: “Of all possible
displacement states a body can assume which satisfy the compatibility equations and
boundary conditions, the state which satisfies the equilibrium equations makes the
potential energy assume a minimum value.”

The potential energy, 17, of an elastic body is given by,

17 = A — S" (2.86)
where A is the strain energy and I” is the work done on the body by external forces.

If the energies are expressed as functions of u, v, and w, then at equilibrium the

principle of the minimum potential energy yields,

H(u,v,w)=A(u,v,w)—'I’(u,v,w) =0 (2.87)

30

The strain energy of a linear elastic body is defined as

A : [aTado (2.88)
.0

where 0 designates the volume of the body. From the stress-strain relations of Eq. (2.77),
and assuming initial strains are present, Eq. (2.88) yields,
A = 1 [57" [R15 d0 - [5T [R].§, .112 (2.89)
2 .0 .0

The work done by external forces is given by,

'11-.- [BTU-(112+[5TU-dr (2.90)
.0 F

where (:0 is the vector of known body forces, O is the vector of surface forces or

fractions, (7 is the vector of displacements and F is the surface where the forces are

applied. Thus,

17: ]§T[R](a—2a,)dr2— [éTU-do— [OTU-df (2.91)

.0 .0 F
Using the variational method, a simple form of the variation of the displacement
field is assumed within each element and the conditions that will minimize the functional

are derived. In this case the functional is given by 17.

Finite Element Equations

The displacement field within an element is given by,

(7 = =[N]Q(e) (2.92)

§<=

31

where [N] is the matrix of shape functions given by Eq. (2.40), and QM is the vector of

nodal displacement degrees of freedom of element ‘e’ given by Eq. (2.39).

The strain vector of Eq. (2.79), using Eqs. (2.83) and (2.84) can be given by,

 

 

a = [B]Q(e) (2.93)
where,
”d I
— O 0
dx
0 9— 0
3y
0 0 .83—
z
[B]— 2- i [N] (2.94)
dy dx
3 a
_ 0 _
dz dx
0 i i
_ dz dy‘
thus,
—b1 0 O C] d] OT
[)2 0 0 C2 d2 0
b3 0 0 C3 d3 0
b4 0 0 C4 (14 O
0 c, 0 b, 0 d,
[B]T ..i 0 C2 0 0 d2 (2 95)
6V 0 c3 0 b3 0 d3 '
0 C4 0 b4 0 d4
0 0 d] 0 bl Cl
0 0 d2 0 b2 CZ
0 0 d3 0 b3 C3
_0 0 d4 0 [24 C4‘

 

 

32

The stresses then can be expressed as

=lRl(5 -Eo)=lRllBlQ‘e’ -lRléo (2.96)

The potential energy can be expressed in terms of the element potential energies

 

 

as
nel
17 = Z 17“?) (2.97)
e=1
where,
l — T T - A T T _
17(9):; I Q“) [B] lRllBlQ‘e’dQ‘eh [ Q“) [B] [R160 40“?)
gm 0(a)
(2.98)
_ J’Q(e) [N17, —dQ(e)_ JQ(e)T[N]T5dr(e)
0(9) r(€)
Thus,
177";le fWB] [1111312099 (272 “3101115, (10(8)
2 e=lg(e) e=l are)
(2.99)
— [[N]T¢ "—dal") [ [N]TOdF(e)
[2(6) [‘(9)
where,
IQl‘
Q=1Qf , (2.100)
lQm.

33

and,

m = degrees of freedom x no. of nodes

(2.101)

The static equilibrium configuration of the body is found by solving the following

necessary condition for the minimization of the functional, 17,

911:6
3Q
This leads to
IKIQ = P

where [K] is the global stiffness matrix,

1K1=§1[K‘6’]

e=l

and P is the global nodal vector
‘ nel _ nel _ nel ‘
424692689248
e=1 e=l e=l
The element stiffness matrix, [Km] , is given by

[168} 113104113149“)
{2(e)

The element load vector due to thermal strains, 3(8) , is given by

2,0 = I 131005.49“)
0(9)

34

(2.102)

(2.103)

(2.104)

(2.105)

(2.106)

(2.107)

The element load vector due to surface forces, PS“) , is given by

163(8) z I [N]TOdF(e) (2.108)
r19)

The element load vector due to body forces, Pf) , is given by

pf): )‘ [N]T p.10“) (2.109)
0(9)

From the theories presented in Sections 2.3.1 and 2.3.2 it can be readily shown

that,
[166)] = V091 [13]T [R][B] (2.110)
and
[1‘
1
(e)
—(e)_EaV . (e)_ . T.‘
P, _(1_2V) (Tm Tmf) [B] <0) (2.111)
0
10.

 

 

where V(e) is the volume of element ‘6’, T (5,2 is the average of the nodal temperatures of

element ‘2’, and T ref is the reference (strain-free) temperature.
(Px Px

Let the body force vector (F = 0),, , and the surface force vector 5 = py , then

(oz [’2

35

Pb“) =7, ) (2.112)

 

 

and,

_. 14(8)
aka—1Q.“ (2.113)

 

 

AS; is the area of face 1-2-3 of element ‘e’. The fourth, eighth and twelfth entries of the

vector in Eq. (2.113) are zero as they correspond N4 which is zero on face 1-2-3. Eq.

(2.113) takes another three forms Corresponding to faces 2, 3, and 4.

36

2.4.3 Reynolds Lubrication Equation
The Reynolds lubrication equation is derived from the Navier—Stokes equations
under the assumptions that:
i. The fluid is assumed to be Newtonian, with direct proportionality between
shear stress and shearing velocity
ii. The fluid is incompressible
iii. The flow is laminar
iv. Inertia and body terms are assumed to be negligible compared to the
viscous terms
v. Variation of pressure across the film is assumed to be negligibly small so
that P = P(x, y) only
vi. The curvature effects are negligible, that is the height of the lubricant film
is much smaller than the length or the width of the body
vii. The viscosity does not depend on pressure but can be affected by
temperature
viii. The motion of the lubricant fluid normal to the surface can be neglected

compared to the motion parallel to the surface.

These assumptions yield the most commonly encountered form of the Reynolds equation,

3[h38P]+a[h38P]18h ah

— —— — —— =—V— — 2.114
dx 1211 dx dy 12/1 dy 2 dy+dt ( )

This equation describes the pressure field of the flow between two bodies with

impermeable walls as shown in Figure 2.9. The first term on the right hand side accounts

37

for the physical wedge and dh/dt =11 accounts for the squeeze film. h is the separation of

the two bodies and ,u is the viscosity of the fluid.

V is the relative velocity given by
V=Va—Vb (2.115)

For a positive pressure field the right hand side of Eq. (2.1 14) must be less than zero.

 

Y

Figure 2.9: Body and coordinate system for Reynolds lubrication equation

Finite Element Eguations

For the derivation of the finite element equations of the Reynolds lubrication
equation the Galerkin method is used. Consider an arbitrary two-dimensional fluid film
as in Figure 2.10 , governed by Eq. (2.114), and with boundary conditions

13:? on r, (2.116)

38

and

3
Vh-i—VP -fi=q 0n r2 (2.117)
1221

where P is the prescribed pressure on F 2, q is the normal flow across the boundary F 2, 17

is the unit outward normal to F2, and I7 = V} .

 

F1 \
O ‘/
X
Figure 2.10: 2D arbitrary body
Eq. (2.114) can be rewritten as
1.3 1 ah dh
V —VP =—V—+— (2.118)
1221 2 dy dt
From assumption (vi) above, Eq. (2.118) can be rewritten as
3
—"—V2P=1Va—h+§-'l (2.119)

1211 2 dy dt

39

Now multiplying by the weight function and integrating over the domain,

3
I—h—VZPwdQI")= ][— V3é+ahjw dale)
”(9)1221 0(6) 2 dy dt

Using Green’s first identity,

[(VP.Vw+mV2P)drz= [(wVP)-h'dr
.0 F

and Green-Gauss theorem

[3% wd.0= 1]}: aydrz+jham (1.0

then Eq. (2.120) becomes

3
[ —h—VP-Vw do“): [[— Vh-a—Q—Qll w]d.(2(e)
0(8) 12,11 0(8) dy dt

(2.120)

(2.121)

(2.122)

(2.123)

If there is no flow across the boundary F 2, that is q = 0, and only the prescribed

pressure boundary condition is true, then I q a) d F (e) = 0 , thus Eq. (2.123) reduces to

rlt’)

are) y are)

Now by approximating,

"o
11
:o
2

40

3
J'[h_vp.vw]dg(9)=%V ] 118—.0) d.(2(e)— [g—Itzwdflm (2.124)

(2.125)

and

leadsto

0(6)

h = h,- N,- (2.126)

g—f=h=li,-N,- (2.127)
wj=Nj (2.128)
3
] MVNi-VN- d0(e)P(e)=-l—V ] (h,N,)_aNf me)
12'” J 2 file) ay

(2.129)

where 13(6) is the vector of nodal pressures of element ‘e.’ Eq. (2.129) can readily be

written as

where,

and,

 

[K(91]P(e) = 1359+ Fife) (2.130)
3
[K99]: I (W) 5N13N1.5N:9N1 me) (2.131)
12;! dx dx dx dx
[2(9)
- 3N.
Fjeile ]' (h,N,-)—J (112(9) (2.132)
2 (e) dy
.0
Fle)=— ] (6.N.)wdr2(e) (2133)
h l l '
{2(9)

41

Assemblage of the element equations leads to

[K]P=F‘

where the global stiffness matrix [K] is

nel

1K1=§[K<e>]

the global load vector F“ is

and,

“01
II

where

m = no. of nodes

 

 

(2.134)

(2.135)

(2.136)

(2.137)

Now from the theories developed in Section 2.3.1 for a linear triangular element,

Eqs. (2.131) to (2.133) can be expressed as

42

 

 

 

 

 

1 12,3
(6‘
31112122 2
2
3h1h2 2
h; 6
2
[Km]:_A_(e,)__(bib,+6ic,)< 3h1-h3 L ..2?
60x12xp J J 6h1h2h3 1
31122113 2
2 2
3"1’73 2
2
3h2h3 (6
. h? .
2 VA”

F119) 6 hile

1

2 1 1

(e)
—(e)__A
h _ 12 1 2 1h,

2.4.4 Incorporation of Boundary Conditions

(2.138)

(2.139)

(2.140)

Afier assembling the element stiffness matrix and the element load vectors the

system of equations can be written as,

[K]a'6=13

(2.141)

However [K] will be singular hence its inverse does not exist and thus the system

Cannot be solved for 5. In the case of static analysis, the physical meaning of this is that

the body is free to undergo unlimited rigid body motion unless some support constraints

43

are provided to keep the body in equilibrium. Thus the prescribed boundary conditions

for a given problem are applied to Eq. (2.141) before solving for d).

If (D j is prescribed as ¢f , the load vector 13 is modified as

*

gag—K ab. (2.142)

1.].
fori = I, 2, m

The rows and columns of [K] corresponding to (D j are made zero, and the diagonal term

is set to one.

K.. =K.. :0 (2.143)

fori= 1, 2, ...,m
and

K j) =1 (2.144)

Finally the prescribed value of d); is inserted into the load vector,

P]. =¢j (2.145)

2.4.5 Guyan Reduction
The Guyan reduction method [20], also known as static condensation, eliminates
the unwanted degrees of freedom of a system, thus reducing its overall size and

consequently decreasing solution time. Given a system,

_.

[K]-X'=P

(2.146)
nxn nxl nxl

where n is the number of degrees of freedom, decomposing it,

44

 

po qu leL pxl}
=<

 

 

 

 

 

 

 

where

q+p=n

(2.147)

(2.148)

Now if K2 is the vector of unwanted degrees of freedom, then the above system can be

reduced to,
[K11 “K121931191?“ = P]
where the reduced matrix, [KW], is given by,
[Kred] = [K11 — K121931191]
The compliance matrix, [C] is given by,
[C] = [Kred ] -1
Also if )‘(2 , the deformation of the rest of the body, is needed then,
X2 = [4931191] 2?,
where the coordinate transformation matrix, [T] is defined as,

1r1=[-K2;‘1<21]

(2.149)

(2.150)

(2.151)

(2.152)

(2.153)

When the stiffness matrix of a static problem is reduced no information is lost. It

can be seen from Eq. (2.150) that all the elements contribute to the reduced matrix.

45

2.4.6 Conjugate Gradient Method

The conjugate gradient method is an iterative method for the solution of systems
of equations. The method cannot compete with direct methods, like solution by matrix
inversion, for small systems. However its strength is revealed with large, sparse systems.
A sparse system is one where most elements of the square matrix are zeros. For such
matrices only the non-zero elements can be stored thus saving storage space. Typically a
matrix is considered sparse if its density is less than 10%. Figure 2.11 shows a graphical
representation of the population of a stiffness matrix obtained from a piston thermal
analysis. The total size of the matrix is 10342, with 11338 non-zero elements. The density
of this matrix is 1.06%. Huge computational effort would be required to invert it, but a

solution for the system can easily be obtained using the conjugate gradient method.

 

270 ~

 

Rows

540 ~

 

810»

 

o 270 ‘ 540 810 1034
Columns nz = 11338

Figure 2.11: Graphical representation of the population of a sparse matrix.

46

Now consider the following function,
f(f)=%iT[A]i—XTI3 (2.154)

The imposed is to find a vector 56 that minimizes the scalar function f (56). The function
is minimized when its gradient is equal to zero, that is
Vf=[A]7c—13=0 (2.155)
This is equivalent to solving
[A]; = 13 (2.156)
The gradient method achieves minimization by iteration. An initial guess for the vector it
is required, 561 , and each subsequent iteration computes a refined solution
56,-“ = 56,- +41% (2.157)
where i is the iteration number, A is the step length and § is the search direction. Eq.
(2.157) must satisfy Eq. (2.156), that is
[A](f,-+zl,-§,-)=E (2.158)

Now let the residual of each iteration be

i=5-[Al‘i (2.159)
then Eq. (2.158) becomes
41415,- =;,. (2.160)
which readily leads to
517?}

(2.161)

47

The search direction is chosen by intuition to be equal to the residual as it gives the
largest negative change in f (if) ,
§,- 2 —Vf = fj- (2.162)
This is known as the method of steepest descent. It is not a very popular method as it
converges relatively slowly. The faster conjugate gradient method utilizes a modified
search direction given by
4.1 = 4+1 +4-4- (2.163)
The constant ,8 is chosen so that two successive search directions are conjugate to each

other, that is

517-11 [4151 = 0 (2'164)
From Eq. (2.163)
(afimrfflda =0 0-165)
which yields,
.7 _
,B- = _ ”1+1 [A131 (2.166)
I 3fiTlAlr-éi

The conjugate gradient method is a very efficient algorithm thus making it very
popular for solving large systems. This method is used in this work for the solution of the
systems of equations developed by the finite element formulation. Figure 2.12 outlines

the conjugate gradient algorithm.

48

   
       

51+1 = F1+1 + fliEi

 

 

 

 

 

(Ll/4151

$714111

fit":—

 

 

 

Figure 2.12: Conjugate gradient method flowchart

49

CHAPTER 3

ENGINE OPERATING PARAMETERS

3.1 Introduction

This chapter introduces the engine operating parameters used in this work. The
derivation of time step, piston position, piston speed, and piston acceleration is achieved
from the engine geometry and engine operating conditions. The theory developed in this
chapter is from Reference [25].

Figure 3.1 shows the terminology used in this work when referring to different

piston areas: the crown, the ring-pack area, the skirt, the pinhole, and the underside.

Crown

Ring-pack
area

Skirt

 

Underside

Figure 3.1: Piston terminology

50

3.2 Time Step

Given a an engine operating at N revolutions per minute (rpm) then the time step,

At, between successive crank angles is given as,

 

 

i=N rev ]X36Odegxlmln (3.1)
At mm lrev 605
thus,
1
At =— (3.2)

3.3 Angular Speed

Given a an engine operating at N revolutions per minute then the angular speed,

a), is given by

 

(1)=N i"- xfixlm‘“ (3.3)
mm lrev 608
thus,
1tN
30 ( )

3.4 Piston Motion

Given the piston assembly shown in Figure 2.1, the piston position, velocity and
acceleration can be readily derived. r is the crankshaft radius, 1 is the connecting rod
length, d is the piston pin offset; 6 is the engine crank angle, to is the angular speed, and (p
is the angle the connecting rod makes with the piston axis. When 6 = O the piston is at the

top dead centre (TDC) at the beginning of the intake stroke. A four-stroke engine requires

51

720 deg. for a cycle, where a cycle has four stokes: intake, compression, expansion and

exhaust.

d

  

ANTI-TH RUST
..... -pX

I

4
THRUST 1
I

 

[12 — (r Sin 6- 692]“2

rcos 19

Figure 3.2: Piston assembly

Figure 3.3 shows a schematic of a piston confined within the cylinder bore. The
piston position at TDC and bottom dead centre (BDC) can be seen. C is the nominal

piston to cylinder bore clearance.

52

THRUST ANTLTHRUST
TOP

DEAD “S """’X
CENTER \ j ; t ,

CYUNDER
BOTTOM §
I - I

BORE
CENTER

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.3: Piston position extremities: bottom dead center and top dead center

Now let
6 = a) t (3.5)
so that
d6
_ = a) 3.6
dt ( )

The piston position yp is given by

{—2 2 1 r 2 d 2
yp= l -—d +r l-cosl9—— 1—[7) [sin9——] (37)

r r

 

53

Eq. (3.7) forces yp to be zero when the piston is at TDC. It can be easily modified to set yp

to be zero at either mid-stroke or BDC.

The piston speed, vp, is given by

 

dyp d6

v = —.__

P d6 dt

(3.8)
rcos6 sin 157-‘1
= no sin 6 + (2 fl) 2
I‘ll—(m (sine—cg)
The piston acceleration, ap, is given by
_ dvp d6
a __ __ .—
P d6 dt
(3.9)
f )

 

2 . _d 2 2 _- - _d
“02 0086+[Lj3 cos 6(sm6 %) 3/2 +[rjcos 9 sm6(s1n6 /r)

(WW-41 4114211543

1

 

 

\

The angle (p can also be readily derived. From Figure 3.2,

Isin¢=rsin6—d (3.10)
thus
¢=sin-l[-———rsmla—d) (3.11)

54

3.4.1 Sample results
Given the engine geometry and operating conditions of Table 3.1, Figure 3.4 to
Figure 3.7 show typical plots of piston position, piston velocity, piston acceleration and

angle {/1 respectively.

Table 3.1: Engine geometry and operating conditions
Crankshaft radius, r: 90.60 mm

 

Connecting rod length, l: 133 mm

 

 

 

 

 

Piston pin offset, d: -0.4 mm
Engine Speed, N: 1000 rpm
0.1 1 l T l l I f
0.03 —————————— —————— ‘ ..... _i ......
E0-06 """" I “‘—r----i— ----- E ————— 1|. _____ E _____ ' _____ _4
C .
§o.04-—-—-——-——7 ..... ; ..... ..... .u .
o.oz--- ----- e ----- ----- 1 ----------
o i i i ' L 1 1
0 90 180 270 360 450 540 630 720
Crank angle [deg]

Figure 3.4: Piston position

55

Acceleration [ms'2]

Velocity [ms’1]
O

N

I
N

is

800

——————————— ————— ————— ——————
----- ----------- ————— -— ~~~~~ 4
..... ............ .....
o ................ I ..... I ___________________ z
-200L ---------- 1 -------------- ; ————— ————————— 4
400 4: ' E : 5 I

 

___________________________________________

 

5.__—__4_____.. .__-__._ ___ __________________

 

0 90 1130 220 360 450 540 630 720

o 90 180 270 360 450 540 630 720

 

 

Crank angle [deg.]

Figure 3.5: Piston velocity

 

 

 

 

Crank angle [deg.]

Figure 3.6: Piston acceleration

56

 

 

 

 

 

0
n 4 4 4| 4 4 2
_ _ r _ 7
_ H _ _ _
H u H n H . 0
T Ilq-lll. lllll fi lllll .llllfilllJllllWl 113
. A _ _ _ 6
_ _ _ _ _ _
_ . _ . r .
_ 7 _ _ _ _
_ m _ _ _ _
_ _ _ _ .
_ _ _ _ . , OJ
TI 1.4 |||||||||||||||||||||||||||||||| 45%
H _ . _ . _
_ _ _ _ . 4m.
_ _ _ . _ _
. _ _ u _ . 0b
TENS-” ........... _ zwm.
“ u H _ _ a
_ _ . _ _ M
1111_-11L 111111111 .1111_11111_1111_1- 11mm
. . _ _ _
. _ _ . _ 2C
. _ _ . r
U u H _ 0
lllllllll .1111. 1. (1.1111111114111118
. _ _ _ _ 1
_ _ . . _ _
_ _ . . _
. _ _ _ . _ 9
_ _ _ _ . . _
_ _ u u H u
. _ _ _ .
.l _ _ ~11. » L 0
4 3. 2 1. 0 1 2. 3. 4.
o o o o mu 0 6.6 m.

Figure 3.7: Angle (p

57

CHAPTER 4
PARAMETERIZATION AND FEA APPROACH FOR THE ASSESSMENT OF

PISTON CHARACTERISTICS

4.1 Introduction

Elastohydrodynamic lubrication, piston dynamics and friction are important
characteristics determining the performance and efficiency of an internal combustion
engine. This chapter presents a finite element analysis on a production piston of a
gasoline engine performed using commercial software, the COSMOSDesignSTAR, and a
comprehensive cylinder-kit simulation software, CASE, to demonstrate the advantages of
using a reduced, parameterized model analysis in the assessment of piston design
characteristics. The fill] piston model is parameterized according to the CASE
specifications. The two are analyzed and compared in the COSMOSDesignSTAR,
considering thermal and mechanical loads. The region of interest is the skirt area on the
thrust and anti-thrust sides of the piston. The results are compared with the CASE results
and a discussion follows on how the piston model simplification approach can help
significantly reduce computation time without losing valuable information on piston
performance. The piston characteristics are evaluated with the reduced model, thus
allowing for a faster analysis and optimization than would be possible with a complete
FEA piston model. The results demonstrate that the simplified model agrees qualitatively
and quantitatively with the full model; the maximum temperature and calculated
distortions on the skirt face are comparable. Computation time difference between the

full model and simplified model is not significant for a single analysis; however it

58

becomes very significant when a cyclical analysis is performed at every crank angle

degree.

4.2 Piston model FE analysis

The finite element analysis is performed in COSMOSDesignSTAR, for both firll
and parameterized piston models. The generated mesh uses tetrahedral elements; Figure
4.1 shows mesh for the full model. A finite element analysis is performed in the CASE as
well. It generates a mesh for the parameterized piston using hexahedral and pentahedral

elements (Figure 4.2). The parameterization of the piston model is shown in Figure 4.3

(a) (b)

Figure 4.1: Full piston model mesh in COSMOSDesignSTAR, (a) top
side, (b) underside

(a) (b)

Figure 4.2: Parameterized piston model mesh in CASE, (a) top side, (b)
underside

and Table 4.1.

 

59

 

 

 

  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

L244

SPA——

  

 

 

 

SW

 

 

Figure 4.3: Parameterized piston dimensions

Table 4.1: Parameterized piston model feature labels

 

PDUk
PHT:
SHT:
CTH:
xpnm
YPDt
XCG:
YCG:
SHTB:
S“k
S“H:
SWQ:
BTH:
PBTH:

PINTHD:

SIR:
SIRC:

piston diameter

piston height

skirt height

mownflmwmmm

x-coordinate of pin location
y-coordinate of pin location
x-coordinate of center of gravity
y-coordinate of center of gravity
tab height

skirt width

lower tab width

upper tab width

boss thickness

pin boss thickness

pin boss width

inner skirt radius

inner skirt radius center

 

60

4.2.1 Thermal Loading

The engine modeled in CASE is a motored optically accessible engine used for
experimental oil film thickness measurements. The engine uses a production piston with
a sapphire cylinder wall.

The engine was first modeled in Ricardo WAVE to obtain the heat transfer above
the piston crown. Figure 4.4 shows pressure traces for the engine, both experimental and
from the WAVE model. The engine was run at 1500 rpm with 2 psi boost. The two
pressure traces have strong agreement except where the exhaust valve opens. The most
crucial part however is at the compression and power strokes. Thus the results from the
WAVE model can be considered reliable. Figure 4.5 shows the temperature of the air in
the cylinder over one cycle. Figure 4.6 shows the heat transfer coefficient above the

piston crown. The heat transfer coefficient is obtained using an advanced combustion

module, the IRIS Model in the WAVE.

 

30 l l I I I T I
----- AVE model
25 , 1 2 I — Experimental ,4

 

 

 

 

N
0

Pressure [bar]
a:

10

 

 

 

 

 

 

0 90 180 270 360 450 540 630 720
Crank angle [deg.]

Figure 4.4: Pressure trace of optical engine

61

Heat transfer coeff.

800 a

 

Temperature [K]

 

 

 

 

 

 

700 ‘ .
600 —————— - - — — — ..... .....
500 ..... L-
400
300 L L i r 4 ' .
0 90 180 270 360 450 540 630 720
Crank angle [deg]

Figure 4.5: In-cylinder air temperature for optical engine

900 r . . , , 4 r
' 1 —- Annand
' '''' Woschni

 

 

 

g 600

400
300
200
100

 

 

 

 

O 90 180 270 360 450 540 630 720
Crank angle [deg.]

Figure 4.6: Heat transfer coefficient above piston crown using
Annand and Woschni correlations

62

WAVE calculates the heat transfer coefficient using both Annand’s [1] and
Woschni’s [55] correlations. Averaging over the cycle, Eq. (4.1), Annand’s correlation
suggests a heat transfer coefficient of 230 W/mzK. Similarly Woschni’s correlation
suggests a value of 170 W/mzK. The average of the two correlations is chosen to be used
as the heat transfer coefficient above the piston crown. The mean temperature, Eq. (4.2),

is found to be 470 K.

— l 47:

8=Zilo hgda (4.1)
.—_ 1 41:
TVWJ‘O hngdB (4.2)

Woschni [56] solved for the heat transfer coefficients governing a piston using
numerical method and the electrolytic tank analog. His results suggest an average heat
transfer coefficient of 460 W/mzK above the crown, in the range of 1000 — 1900 W/mzK
in the ring groove area, 220 — 930 W/mZK at the land area, 550 — 720 W/mZK at the
piston skirt and 100 — 1300 W/mzK on the piston underside. Similar values are used by Li
[35].

In the study of the parameterized piston, the heat transfer coefficients are assumed
constant in each area and within Woschni’s range (Figure 4.7) and having values as in

Table 4.2.

63

 

 

 

 

'1L__T__1_r

 

 

 

 

 

Table 4.2: Heat transfer coefficient and ambient temperature values

 

12/f

11/ \

Figure 4.7: Heat transfer coefficients for parameterized piston

 

 

for parameterized piston
g = 200 W/mzK T, = 470 K
h,=1000W/m2K T,=314K
h; = 700 W/mZK r2 = 310 K
h = 250 W/mZK T3 = 303 K

 

 

The same values are used for the full model analysis with the addition of a heat

transfer coefficient at the land area as well as at the pinhole of 400 W/mzK and 800

W/mzK respectively.

The thermal analysis was performed using three different sets of heat transfer
coefficients: the ones of Table 4.2, an increase by 15% and a decrease them by 15%. This
was done to assess the sensitivity of the piston to changes in the heat transfer coefficients.

The full piston model temperature distributions for the three cases are shown in Figure

4- 8, Figure 4.9, and Figure 4.10 respectively.

64

3 247€+U|32
"I ‘ '3 2355-4102
2236'002

3
3 2’1 28+DCIZ

    
  

3 200mm:
3 188124002
3 1778‘002
3 1555:4002
3 1536+002
3 142E+DOZ

3 1 3069002

Figure 4.8: Full piston model temperature distribution with boundary conditions of
Table 4.2 (temperature in Kelvin)

Temp
3 2708+002

' 253644302

  

3247844302
,, _3 235e+002
- , 32232o002
' . 321213.002
7 3200e+002
3188:3002
3177134002
3155124002

: 3 15313‘002

 

3 142891302

3 130942102

 

Figure 4.9: Full piston model temperature distribution with heat transfer coefficient
values increased by 15% (temperature in Kelvin)

Temp
3 2709002
3 2586+002

3 247E602

    

 

132¥em02
' ‘:.3223e0002
3: 32129»002
320069002
3188e+002
3 'l 77e+002
' 316$e+002
3153e‘UU2

“ 3142e+002

   

3 130134-002

Figure 4.10: Full piston model temperature distribution with heat transfer coefiicient
values decreased by 15% (temperature in Kelvin)

Performing the thermal analysis on the parameterized piston model, it can be seen
that the temperature distribution follows a very similar pattern as the full piston model

(Figure 4.11, Figure 4.12, and Figure 4.13).

Temp
3 30091002

3 28883002

3 2'1 39002
3 20139002
3 1888+UU2
3 1758"]:12

3 16 393002

 

3 lSUe+002

Figure 4.11: Parameterized piston model temperature distribution with boundary
conditions of Table 4.2 (temperature in Kelvin)

66

Lo
L_'.
1_.
'3
'1
$
'2‘
1:.
r J

   

3 2756.002
3 2636+002
, 3 2509002
3 230e+002
’ 32259002
é 321324002
3 200.:«002
318Eie+EIIJ2
‘ , 31754002

' 3153e+002
3 1509002

Figure 4.12: Parameterized piston model temperature distribution with heat transfer
coefficient values increased by 15% (temperature in Kelvin)

Temp
3 3IJUe+DO2
3 2889+1302

3 2.7581002

 

" . 3 26319002
3 2504002
3 23013002
3 225e+002
3 21364002
3 2DOe+UO2
1. 3 113013.002
‘ 3 1759002
I 3 153e+002

 

3 TSDe+UOZ

 

Figure 4.13: Parameterized piston model temperature distribution with heat transfer
coefficient values decreased by 15% (temperature in Kelvin)

Table 4.3 summarizes the maximum and minimum temperatures for each of the
above cases. The trends are the same for both models. The parameterized piston model,
though, has less surface area exposed to convection than the full piston model.

Consequently it has a maximum temperature of about 4 K higher than the full model.

67

Also the minimum temperature of the parameterized model is about 2 K higher. These
results demonstrate that the temperature distribution on the piston is not greatly affected
by the heat transfer coefficient values. The temperature difference of nodal values is in
the range of 2 — 4 K which can be acceptable since the complexity of piston thermal

loading is simplified as in Figure 4.7 and Table 4.2. Thus the assumed values of Table 4.2

can be considered valid.

Table 4.3: Maximum and minimum temperature for the
three sets of heat transfer coefficients

 

 

 

h Tmax (K) Tmin (K) AT (K)
Full Piston Model
Table 4.2 325.5 314.2 11.3
+15 % 326.2 313.6 12.6
~15 °/o 324.7 314.7 10.0

 

 

 

 

Parameterized Piston model

 

 

Table 4.2 329.2 315.9 13.3
+15 % 330.1 315.3 14.8
-15 % 328.3 316.6 11.7

 

 

 

4.2.2 Thermal Loading Deformation

The temperature distribution on the piston causes thermal strains which
consequently deform the body. Considering only deformation due to thermal strains
(Figure 4.14 and Figure 4.15), the two models are not in strong agreement quantitatively;
but qualitatively they both deform outwards, away from the center axis. The fill] model
has a maximum deformation of 102 um at the tab edge and the parameterized 78 pm at

the same location. However, considering the entire skirt face, the average deformation is

68

58 um for the full model and 57 pm for the parameterized. Thus on a mean basis the two

models have a strong agreement.

LlFTES
1 03064301

9 4429-002

  

13 503.2002
‘ g, r 7 725.2002
0 857.2002
5 0089.002
5 150e.002
4 2926—002
3 433e.002
2 575e002
1 7173002
8 50312003

   

0 00020000

Figure 4.14: Full piston model deformation due to thermal loading (URES: resultant
deformation in mm, deformation scale 1:50)

URES
'1 03013001
9 4428-002

  

8 58313-002

 

"1 L’ 1 2 22512002
13 00719-002
: 600313-002
515019002
4 29212002
3 3312.002
2 57512002
1 71712002
8 583e-003

 

0 UDDe+CIIJU

Figure 4.15: Parameterized piston model deformation due to thermal loading (URES:
resultant deformation in mm, deformation scale 1:50)

4.2.3 Mechanical Loading

The piston experiences mechanical loads that force it to deform. In the analysis
three loads are considered as well as a constraint. Ambient pressure (100 kPa) is applied
above the crown to simulate the effects of combustion gases, and on the piston skirt to

2 is applied to

simulate the effects of lubrication pressure. A body load of 9.81 ms'
account for inertial effects. The pinhole is restrained (Figure 4.16). The same loads are

applied to the parameterized piston model.

 

Figure 4.16: Mechanical loading, red arrows: atmospheric pressure, green arrows:
constraint, body load not shown

Figure 4.17 and Figure 4.18 show the deformation due to mechanical loading for
both the fill] piston model and the parameterized model respectively. The two
demonstrate exactly the same behavior. The full model has a maximum deformation of
5.98 um, and the parameterized model deforms 6.04 um towards the piston center axis.

The two models have a strong agreement both quantitatively and qualitatively.

7O

5 03 33003
.7, i 4 5756-003
,_- I 4 087134303
3 5586-003
3 05135.13le
2 5425.003
2 0 3 312-003-
1 5255003 1
I 0176-1303

5 136' 3134304

 

U |:l[|l:le+1:ll:ll]

Figure 4.17: Full piston model deformation due to mechanical loading
(URES: resultant deformation in mm, deformation scale 1:50)

LlRES
5 101191103
S 59213—003

S 08'3e~0[I3

 

' 4 57564303

=10:

' 4 [16780133

    

3 0506-003
2 5426-003
2 0338-003
~« 1 52664303
1 0178-003

' ,5 5003mm

 

0 C'UU‘?+000

Figure 4.18: Parameterized piston model deformation due to mechanical loading
(URES: resultant deformation in mm, deformation scale 1:50)

4.2.4 Combined Thermal and Mechanical Loading
This section examines the effects on piston skirt deformation of both thermal and

mechanical loads applied together. Since the lubrication pressure on the skirt is not

always uniform, influenced by piston tilting, piston lateral motion and bore distortion,

different skirt load cases are considered in examining the overall skirt behavior.

0 Case 1: Uniform normal pressure on skirt (Figure 4.19a)

0 Case 2: Directional uniform pressure on skirt (Figure 4.1%). Pressure is in
the z-direction

0 Case 3: Stepped uniform normal pressure on skirt (Figure 4.19c).
Maximum pressure is at the top of the thrust side.

0 Case 4: Stepped uniform normal pressure on skirt (Figure 4.19d).
Minimum pressure is at the top of the thrust side.

0 Case 5: Stepped uniform normal pressure on skirt (Figure 4.19c).
Maximum pressure is at the center of the skirt.

0 Case 6: Stepped uniform normal pressure on skirt (Figure 4.191“).
Maximum pressure is at the right half part of the piston, looking at it from
the thrust side.

The rest of the loads and constraints are as in Figure 4.16.

In Figure 4.20 to Figure 4.31, the deformation results are shown for each of the

above load cases. The thrust side is shown. The resultant deformation (URES) is in mm.

The deformation scale is 1:50.

72

      
 

Thrust Sid.

 

i N Thrust ha.

 

 

©,-——> :1E+5 Pa —->:1E+4Pa Qu—D :1E+3 P-
W

-—> :1E+5PI

 

 

 

Figure 4.19: Pressure loads on skirt

URE'S
9 3008-002
8 9336-002

8 1573—002

r 17 3506-002

6 5336-002

1

5 7176002
4 90025002
4 08313-002
3 257E«UU'2
2 4509-002
1 633184302

8 IEII’EuOUB

 

0 000690013

Figure 4.20: Full piston model deformation due to load case 1

73

LIFJES
0 3002.002
3 90312002
3 1575.002
1 ii- _ 7 350122002
' ‘ a 533e.002
' f 5 71 715.002
4 900.2002
4 083e-002
3 20712002
2 4505.002
1 6339.002

‘ 51 15712-003

   

Figure 4.21: Parameterized piston model deformation due to load case 1

URES

9 sane—002

8 133712002

8 0332-002

_ 7 230e-002
15 42712-002

  

5 5238-002
4 820(2-002
4 0176-002
3 2136-002
2 41 DIE-002
1 5078-002
8 0336-003
0 UIZIIJe+000

 

Figure 4.22: Full piston model deformation due to load case 2

URES

Q 640e—002
8 8378-002
8 0336-002

4 1 7 23013-002

‘ 6 427.2002
5 62312-002
4 8208-002
4 0178-002
3 21 32-002
2 41 file-1302
1 son—002
8 0338-003
0 00Cle+000

Figure 4.23: Parameterized piston model deformation due to load case 2

    

LIRES

9 4006-002
8 81 78-002
7 8338—002

62678-002
5 4836-002
4 7008-002
3 9178-002
3 133e-002
2 3508-002
1 5579-002
7 8338-003
0 0009000

 

Figure 4.24: Full piston model deformation due to load case 3

75

URES
9 40013-01112
8 61 79-0132

  

7 8339-002
, ‘. 7 050.3002
15 26713002
1 Y 5 4836-002
4 70012002
3 917.2002
3 133-2002
2 35012002
1 56713-002
7 033.2003

 

0 0006+000

 

Figure 4.25: Parameterized piston model deformation due to load case 3

LIFE?)

9 620e-002

8 3186-002

1.‘ 8 0178-002
1. 7 21513-111132
' 5 41313002
5 6128-002

4 8108-002

4 Ulih‘Be-OIJZ‘

3 2078-002

2' 4056-002

1 5036-002

8 0176-003

0 000239000

Figure 4.26: Full piston model deformation due to load case 4

0053
9 5205.002
8 81813-1102
8 0179.002
_ 7 2155.002

‘ 5- 4135-002
5 151212.002
" ' 4 0105.002
{' 4 [1083002
3 2075002
2 4052.002
. , 1 5035.002
‘ ‘ 8 0175003

 

0 0006 #000

 

Figure 4.27: Parameterized piston model deformation due to load case 4

URES

9.600e-002
8 3008-002
0 0005002
F1117 20013—002

‘ -; 6 40015002

  

5 6005002
4 60013-0132
4 0002-002
3 2006—002
2 400e-002
1 60013-002

8 00013-003

 

0 000e9000

 

Figure 4.28: Full piston model deformation due to load case 5

LIRES
9 6006—002
8 3009 -002

8 0008—002

 

7 2005002
13 4005002
5 0005002
4 0005002
4 0005002
3.2005002
2 4005002
1 5005002
0 0005003

0 00013-1000

 

Figure 4.29: Parameterized piston model deformation due to load case 5

UF’ES

9 5006-002
8 7088-002
7 9178-002

 

7 1 258-002

11.1..

6 3336-002
5 5428-002
4 7506-002
3 9588-002
316713-002
2 3758-002
1 5836-002
7 917e-003
0 0006+000

 

Figure 4.30: Full piston model deformation due to load case 6

78

URES
9 50012-002
S 7086-0132
" 7 9175002
' I 1 'L 7 125-94302
‘ 1’ 1 5 3336-002
‘ i 5 5429-002
4 7505002
3 9588-002
3 1676-002
‘ 2 3758-002
‘ 1 5836-002
7 5175003

 

0 000¢+000

 

Figure 4.31: Parameterized piston model deformation due to load case 6

Considering Figure 4.21 to Figure 4.31 and Table 4.4, it can be concluded that the
deformation on both the firll model and the parameterized model follows the same trends.
The two models have a qualitative agreement for all six load cases. The maximum
defamation occurs around the tab area. The thermal strains, as discussed earlier, lead to a
poor quantitative agreement, about 2 pm for each load case. The mean deformation has a
difference of about 1 um. The maximum deformation occurs at the tab edge of the anti-
thrust side for both models. It occurs on this side because the pin offset is towards the
thrust side. Consequently the restrained nodes, at the pinhole, are further from the anti-

thrust side, which allows it higher deformation.

79

Table 4.4: Maximum and mean deformations on skirt

 

 

 

 

 

Parameterized
Skirt Full Model Model
Load Deformation (um)

Max Mean Max Mean

Case 1 97.0 53.0 73.8 51.9
Case 2 96.4 52.9 72.9 51.7
Case 3 94.0 51.7 76.2 52.6
Case 4 96.1 51.4 77.2 52.3
Case 5 95.6 51.6 76.9 52.6
Case 6 94.6 51.4 75.9 52.4

 

 

 

 

4.2.5 Computation Time

An independent finite element solver running under MATLAB was used to

compare computation time between the full model and the parameterized model (Table

4.5). The mesh (4-node tetrahedral elements) was imported from COSMOS for both

models. The analysis was run on a Pentium 4 3.0 GHz.

Table 4.5: Mesh and computation time comparison

 

 

 

 

 

Parameter Full Parameterized
Model Model

No. of Elements 25270 6590

No. of Nodes 6339 1959
Thermal Loading

CPU Time (sec.) 30'” 3'76
Mechanical

Loading 489.55 17.03

CPU Time (sec.)

 

 

 

In the case of the thermal loading analysis the solution time for the full model is

greater by a factor of 8, compared to the parameterized model. Similarly in the case of the

mechanical loading analysis the full model solution is slower by a factor of 29.

Consequently running a full model analysis for one complete cycle would require 97.9

80

hours compared to 3.4 hours for the parameterized model. Thus the parameterized model
is more cost-effective. Also it allows for a faster optimization of the design
characteristics.

The mesh generated in CASE (pentahedral and hexahedral elements) for the
parameterized model has 1620 elements with 6169 nodes. CASE uses a FORTRAN
based finite element solver which required approximately 20 seconds for the solution of

the thermal loading and the mechanical loading analyses.

4.3 CASE Analysis
4.3.1 Finite Element Analysis

Figure 4.32 shows the FEA results obtained from CASE for thermal loading. The
temperature distribution ranges from 318 K to 337 K. It is slightly higher than the
COSMOS result; the CASE uses a coarser mesh, also the solvers used are different. The

COSMOS uses iterative solvers, whereas the CASE uses Gauss elimination.

» ‘ if“... . »
F‘NW‘VL 202:“-

 

Figure 4.32: Parameterized piston model temperature distribution from CASE
(temperature in degrees Celsius)

81

Figure 4.33 shows the FEA results obtained from CASE for the combined
mechanical and thermal loading. These results, though, do not include the ambient
pressure loading on the piston skirt. CASE adds pressure loads on the skirt during the
cyclic analysis when there is a lubrication pressure. The deformation distribution follows

the same pattern as the COSMOS results with a maximum deformation of 125 um.

 

Figure 4.33: Parameterized piston model deformation obtained from CASE
(deformation in mm)

4.3.2 Piston Characteristics
This section demonstrates some of the capabilities of CASE. Only some of the

results are shown here: oil film thickness, friction force and piston dynamics.

Oil Film Thickness
Figure 4.34 and Figure 4.35 show the oil film thickness distribution on the skirt at
33 deg. on the thrust and anti-thrust sides respectively. From these, the piston orientation

can be deduced. At 33 deg., the top part of the piston tilts towards the thrust side and the

82

bottom part towards the anti-thrust. The oil film thickness distribution follows the pattern
of the deformed skirt profile as in Figure 4.20. Also it is in the range of 500 pm, which is
the clearance between the nominal bore diameter and nominal piston diameter. Nominal
clearances are usually in the range of 50-70 pm. Here, the clearance is very large because

of the specifically built sapphire cylinder.

Friction Force

Figure 4.36 shows the friction force between piston and cylinder wall. A positive
value designates contact with cylinder on the anti-thrust side and negative value
designates contact on the thrust side. The four contact points considered are shown in

Figure 4.37. When the friction force is zero, the piston is floating in the cylinder.

Thickness

[um]

    

Distance from bottom of oil ring groove [mm]

-10
Distance from skirt center axis [mm]

Figure 4.34: Skirt oil film thickness on thrust side at 33 CAD

83

Thickness

[11m]

   

450

 

‘ 400
; 350
, 5- 300
250
;2oo

 

 

  

Distance from bottom of oil ring groove [mm]

   

-270;7 ’ -1o 0 1o 7
Distance from skirt center axis [mm]

Figure 4.35: Skirt oil film thickness on anti-thrust side at 33 CAD

800

 

 

 

 

 

, Topcontactpoint , .

600 ..... 3 . B°§8mflfitflpflm ‘
4oo----

-Z- 200 ,, l-

‘” \i I i

9

O . _____ »V 7 ,,., 773

“E o ‘ 1 r‘ . ‘ ‘ :

Q . ,

22-200» - ”1‘ ,

Uh. 1 . I . . ‘
400- . ------ ; - .‘
-6oo~ 1 ~~~~~~ rrrrrr
-800 i i i L l i I

0 90 180 270 360 450 540 630 720

Crank angle [deg.]

Figure 4.36: Friction force between piston and cylinder liner

84

Thrust
Side

Top contact points

/

 

 

\
g

 

 

 

 

#
5%

é

 

 

 

 

 

   

\ Bottom contact point

 

 

Antithrust
Side

Figure 4.37: Contact points of piston with cylinder wall

Piston Dynamics

Figure 4.38 shows the piston rotation over a cycle. At 33 deg. the piston has a

negative rotation, which means that the top part of the piston tilts towards the thrust side.

This agrees with the oil film thickness distribution.

Figure 4.39 shows the piston pin lateral motion. A positive value designates
movement of the piston pin towards the anti-thrust side. Negative values show movement
towards the thrust side. Piston pin lateral motion affects contact forces. Considering

Figure 4.36 and Figure 4.39, it can be seen that when the pin lateral position is at the

maximum the friction force peaks.

85

 

 

 

 

 

 

 

 

 

T1111

 

 

   

 

180 270 360 450 540 630 720

90

 

 

004
0.031-———-——-_- _ ,

m
0.
Gm: 5:32 .555

123
0.0.0.
0.n.un_u

Crank angle [deg]

Cyclic piston rotation about pin

Figure 4.38:

 

 

 

 

 

 

 

 

 

 

 

 

 

2 4
o. a. 0.
FE: coEmoa Ed

180 270 360 450 540 630 720

90

Crank angle [deg.]

Figure 4.39: Cyclic piston pin lateral motion

86

4.4 Experimental Results

The experimental results1 for oil film thickness from the motored sapphire bore
engine are inconclusive. Very few results are available and although they demonstrate a
qualitative agreement with the simulated results they do not agree quantitatively. The
experimental results are obtained using the laser induced fluorescence method. Figure
4.40 shows the experimental oil film thickness measurement at the TDC of the intake
stroke. Figure 4.41 shows the predicted oil film thickness by CASE at the same stroke
and location. Both experimental and simulated results show a thicker oil film near the oil
ring groove, which gets thinner away from it and towards the piston central axis. The
quantitative discrepancy can be attributed to three factors.

At this preliminary stage the thermoelastic and lubrication behavior of the
sapphire cylinder is still not well understood to be applied accurately to the piston model.
In the CASE simulation a fiilly flooded elastohydrodynamic model is assumed, however
experimental results show some irregularities which indicate the potential existence of
the partially flooded condition. Finally in the specific simulation the piston-to-bore
clearance is independent of the elastohydrodynamic behavior of the piston. The piston-to-
bore clearance is calculated from the cold-stage piston skirt profile rather than the crank
angle dependent deformed skirt profile which leads to bigger clearances. Consequently
these factors are under consideration and together with the acquisition of more

experimental results will be presented at a later stage.

 

 

1 Experimental oil film thickness measurements were provided by Andrew Fedewa of Mid-Michigan
Research sponsored by the US Army.

87

Thickness
um

 

 

 

Distance from top of piston [mm]

  

-6 -5 -4 -3 -2 -1 0
Distance from skirt centre axis [mm]

Figure 4.40: Experimental oil film thickness measurement at TDC of intake stroke

Thickness

 

Distance from top of piston [mm]

 

-6 -5 -4 -3 -2 -1 0 1
Distance from skirt center axis [mm]

Figure 4.41: Simulated oil film thickness prediction at TDC of intake stroke

CHAPTER 5

PISTON FINITE ELEMENT ANALYSIS: THE FULL MODEL APPROACH

5.1 Introduction

As was shown in the last chapter, a parameterized piston model can be used as
guidance in assessing the piston characteristics. Once that guide line is obtained, what
can be done with it? The answer is proposed in this chapter. The goal is to develop a
piston program that will supplement the existing one in CASE. The new program will be
able to perform an analysis for a full piston model. The optimal design characteristics for
the piston obtained from the parameterized model analysis would be used to modify the
computer aided design (CAD) model of the full piston. The optimized full model then
would be analyzed to ensure for optimal piston characteristics.

The first module of the new piston analysis program is introduced here. This
module performs the finite element analysis of the fill piston over the cycle, thus called
piston finite element analysis, from hereon referred to as n-fea. The Greek letter ‘7r’

appears from the Greek translation of piston, ‘marévz’.

5.2 Preprocessing

Unlike the existing PISTON program in CASE, n-fea is not a standalone

program; some preprocessing is required for the full piston model.

89

5.2.1 CAD Model

The current version of 1r-fea requires the piston model to be oriented as shown in

Figure 5.1.

 

Figure 5.1: Piston orientation

The x-axis should pass through the skirt face with the anti-thrust side being in the
positive y-z plane. The y-axis should be parallel to the piston axis and the z-axis should be
parallel to the pinhole axis. The origin of the piston’s global coordinate system should be
at the centre of the crown. The orientations of the piston can be very easily changed in

any CAD sofiware. For this work SolidWorks is the software of choice.

90

5.2.2 Mesh and boundary conditions

The next step is to mesh the model and apply the boundary conditions for both the
thermal and static analysis. The current version of u-fea supports mesh data files created
in COSMOSDesignSTAR.

COSMOSDesignSTAR uses only linear and quadratic tetrahedral elements for
meshing. The linear tetrahedral elements do not have the accuracy of the quadratic
tetrahedral elements as they do not conform to the geometry as well due to their linearity.
Quadratic elements however, lead to very large systems of equations. Thus considering
limitations in computer memory, the linear tetrahedral elements were chosen for this
study. n—fea can be easily modified to support different type of elements.

The boundary conditions and restraints are applied in COSMOSDesignSTAR as
well. Rather than applying the numeric value of the boundary condition in
COSMOSDesignSTAR, the piston faces affected by a specific boundary are flagged.
This allows for n-fea to recognize faces and apply the right boundary condition. Also the
user can run different simulations with different boundary conditions without having to
go back to COSMOSDesignSTAR.

Table 5.1 shows the flags used for the thermal analysis. A convective boundary is

applied to the faces. The heat transfer coefficient value is assigned as a flag for each area.

Table 5.1: Flags for thermal analysis

 

 

 

Piston area COSMOSDesignS TAR Flag
Crown h = 1 Wm'ZK'l
Ring-pack area h = 2 Wm'ZK'l
Skirt area (thrust and anti-thrust sides) h = 3 Wm'ZK'l
Underside h = 4 Wm'zK'l

 

91

Table 5.2 shows the flags used for the static analysis. A uniform pressure flag is
assigned to the crown, thrust and anti-thrust sides faces. A fixed restraint is applied to the

pinhole and a body force is the positive y-direction is applied to the whole piston.

Table 5.2: Flags for static analysis

 

 

Piston area COSMOSDesignS TAR Flag
Crown P = 1 Pa
Skirt area (thrust side) P = 2 Pa
Skirt area (anti-thrust side) P = 3 Pa
Pinhole Restraint = Fixed
Body Body force: y-dir = 9.81 ms'2

 

 

COSMOSDesignSTAR outputs two data files, one for the thermal analysis and
one for the static analysis. Each data file contains the nodal coordinates, the connectivity

matrix, and the boundary conditions. The format of the data files is described in reference

[8].

5.3 Piston Finite Element Analysis (n—fea) program

After the piston model is preprocessed the 1r-fea program takes over to perform
the cyclic analysis. Figure 5.2 shows an overview of the 1r—fea program flow. Detailed
flowcharts of the n-fea program and of its subroutines are shown in Appendix A.

The n-fea program starts by loading the mesh data files exported from
COSMOSDesignSTAR. Then it loads the file that contains the piston material properties,
engine geometry, boundary conditions and pressure trace. After that it calculates the
piston motion as described in Section 3.4. Following this it calls the FEMmeshread
subroutine. Here the mesh data files are read. This subroutine returns the nodal

coordinates, connectivity matrix, and arrays containing the boundary conditions with

92

their associated faces. Then it calls the FEMthermal subroutine. Here the thermal finite

element analysis is performed as described in Section 2.4.1.

Geometry and orientation in SolidWorks
Mesh and flag faces in COSMOSDesignStar

Calculate piston temperature distribution

Calculate piston deformation due to
thermal loading,
unit body force loading,
and unit combustion gas pressure loading

Calculate oil film thickness, skirt lubrication
and contact pressures, forces, and moments

Calculate piston deformation due to skirt
pressure loading
and piston total deformation

Calculate strains: ex, y 82, My 7x2 7y: seq

Calculate stresses: ax. 0,, oz, Txy sz' ryz,

01' 02' 03' ”von' 0'int

 

Figure 5.2: Overview of the n-fea program flow

93

FEMthermal returns an array containing the nodal temperatures, as well as an array
containing the element characteristics, that is, volume and shape function derivatives.
This element characteristic array is used by the rest of the subroutines that require
element characteristics information in the calculations.

This speeds up the whole program as element characteristics are calculated only
once. Following the FEMthermal subroutine the FEMstatic subroutine is called. Here the
deformations due to the thermal strains, d,, due to a body force, b, = 9.81 ms'z, db“, and
due to a uniform pressure, Pu = 100 kPa, on the crown, dgu, are calculated as described in
Section 2.4.2. Deformation is assumed to vary linearly with body force and pressure
loading thus dbu and dgu are used to calculate the actual deformation at each crank angle
over the cyclic analysis. FEMstatic also returns the connectivity matrix for both the skirt
thrust and anti-thrust sides. After this, u-fea prepares an array with the node numbers of
the nodes on the skirt in ascending order. This array is used as an input for FEMGuyan
subroutine where the skirt compliance matrix is calculated as well as the coordinate
transformation matrix. The method is described in Section 2.4.5. The Guyan reduction
method allows for the reduction of the system just to the nodes on the skirt, thus allowing
for a faster solution. The next step is to transform the skirt surface (bearing area) of the
thrust and anti-thrust sides from the three-dimensional global coordinate system to a two-
dimensional local one (Figure 5.3). This is done under the assumption that the skirt
curvature is negligible compared to the skirt circumferential length, thus allowing for the

solution of the Reynolds equation as described in Section 2.4.3.

94

 

 

 

XI

,tYS

 

Figure 5.3: 2D skirt profile

Now, with all these prepared, the cyclic analysis begins. The time step is set at
one crank angle and is controlled by the external loop. The deformations due to

combustion gas pressure, Eq. (5.1), and piston axial acceleration, Eq. (5.2), are

 

 

calculated.
P (t)
dg(x,y,z,t)= 81:, dgu(x,y,z) (5.1)
a I
db(x,y,z,t)= pb( )dbu (x,y,z) (5.2)

Following this an iterative loop starts, the internal loop, to calculate the total
pressure on the skirt. An initial guess for the total pressure is required at t = 0. To make
the initial guess simple it is suggested to start from zero pressure, P,(t = 0) = 0. Within the
internal loop the FEMstaticLUB subroutine is called, where the load vector due to

pressure loading on the skirt is built. Using the reduced stiffness matrix form FEMGuyan

95

subroutine the deformation on the skirt due to hydrodynamic pressure loading, d,;,, is

calculated. Having all these, the total deformation on the skirt, d,, is calculated.

ds,m (x,y,z,t)=d,,m (x’y’z)+db,m (x’y’z’t)

(5.3)
+dg,m (x,y,z,t)+dsh,m (x,y,z,t)
where,
1 thrust side
m = (5.4)
2 anti —thrust side

Given the piston orientation described in Section 5.2.1 the lubricant film thickness can be

well approximated by,

hm (t)=C+(R—"ds,m(x,y=c0nst,z,t)|l) (5.5)
where C is the nominal clearance between cylinder bore and piston, R is the piston radius,
and II d3”, (x, y = const,z,t) H is the magnitude of the deformation vector evaluated from

x = 0, z = 0 and at constant y. The squeeze film is given by

km(t)=hm(t)-:tm(t—l) (5.6)

 

If the lubricant film thickness, hm, becomes smaller than the skirt waviness, Q,
then wavy contact occurs between the cylinder wall and the skirt. The wavy deformation,

(5,", is given by

.Q—h ' h [2
5m={ ’"’ 'f "‘< (5.7)

0, ifhm>.(2

Zhu et al [58, 59] , using Johnson’s model (1985) for a blunt wedge against a plane,
proposed a formula for the mean solid-to-solid contact pressure between an aluminum

alloy piston and a cast iron cylinder wall:

96

Pam =(5.464x1013)6,,,‘-°552 (5.8)

This relation is used in this work.

Following this the ReEq subroutine is called. Here the Reynolds lubrication
equation is solved for the hydrodynamic pressure, P1,, on the skirt area as described in
Section 2.4.3. P1, = 0 is used as the boundary condition for the Reynolds equation. Once
the hydrodynamic pressure is calculated, the half-Sommerfeld condition, Eq. (5.9), is

applied as conventional lubricants cannot withstand negative pressures and cavitate [32].
P1,... = P1,... -x(P1,m) (5.9)

where 1(Ph, m) is a switch function,

1, ifpmzo

Z(Ph,m)={0, if PM <0 (5.10)

The total pressure, P,, can then be calculated:
Pt,m=Ph,m+Pc,m (5'11)
The program checks for pressure convergence. The pressure residual, Pm, is calculated

by,

P... = Haj—74'
31

2|

where i is the internal loop iteration number. If the pressure residual between successive

(5.12)

 

internal loop iterations goes below a predefined tolerance (usually 103), then the internal

loop exits to the external loop. Otherwise the pressure is modified using underrelaxation:

97

P. =PX“ +wr(P."—I:"‘) (513)

where a), is the underrelaxation factor. a), begins at 0.5 and if successive pressure
residuals increase, then a), is decreased by 0.1 up to 0.1. P), and PC are very sensitive to
the deformations thus numerical stability becomes more difficult to achieve as the
lubricant film thickness decreases.

Once the pressure is converged, the hydrodynamic, contact and friction forces and

moment acting on the piston can be calculated (Figure 5.4).

THRUST ANTI-THRUST

L
.4

J
l

Ki M lFf
4

 

 

.. .._..-.—._

 

 

 

 

 

 

 

 

 

 

 

Figure 5.4: Forces and moments acting on piston

The hydrodynamic force is given by integrating the hydrodynamic pressure over

the two-dimensional skirt area (Figure 5.3),

Fh,m = H Ph,m (xs,m’ys,m)dxs,mdys,m (5'14)
A

s,m

Similarly the hydrodynamic moment about the wrist pin is given by,

98

Mh,m : H Ph,m (xs,m’ys,m)'(a_ys,m)dxs,mdys,m (5'15)
A

The total hydrodynamic force F), is given by,
Fh :Fh,l_F/1,2 (5.16)
and the total hydrodynamic moment, M1,, is given by,

Mthh,1‘Mh,2 (5-17)

The shear stress acting on the skirt can be expressed as,

r =———P-+—"—’-——’— (5.18)

The hydrodynamic friction and its moment about the wrist can be calculated as follows:

th,m = fl Tm (xs,m’ys,m)dxs,mdys,m (5'19)
AS m

th,m = “- Tm (xs,m’ys,m)'(xs,m 'Cp )dxs,mdys,m (5'20)
A5 I"

The total hydrodynamic friction, F,;,, and total hydrodynamic moment, M,;,, are expressed

as follows

th =th’] +th’2 (5.21)
th = Mflhl _Mflt,2 (5.22)

In a similar way the total contact force, F6, and total contact moment, MC, are

calculated:

99

:49; Pcm x3 m’yS m)dxs,mdys,m (5.23)

M m = H Pc,m (xs,m2ys,m)'(a—ys,m )dxs,mdys,m (524)
PC = FCJ _Fc,2 (5.25)
Me = Mc,1—Mc,2 (5.26)

Given a fiiction coefficient ,uf, between the cylinder liner and the piston skirt, then the

contact friction force, F k, and contact friction moment, Mfc, are expressed as follows:

ch,m =—I_v_ Pl Iflch, m(x x,s m’ys ,m)dxs,deS,m (5'27)
vp AS, 771
Mfc,m :I—': :I AIflch ,m(xs ,m’ys ,m) (xs,m -Cp)dxs,mdys,m (5°28)
ch = ch,l + Ffflz (5.29)
Mfc = Mfc,l _Mfc,2 (5.30)

Now the total normal force on the piston due to hydrodynamic and contact

pressures is given by,

EthJch (5.31)
the total friction is given by,
Ff=th+FfC (5.32)
and the total moment, M, is given by,
M=Mh+MC+Mjh+MfC (5.33)

100

Following this the total deformation of the piston is calculated. Using the
coordinate transformation matrix (Section 2.4.5), [T], the deformation for the whole
piston, db, due to the pressure loading on the skirt is calculated. The total piston

deformation, d, is given by,

d(x,y,z,t) = d, (x,y,z)+db (x,y,z,t)

(5.34)
+dg (x,y,z,t)+dh (x,y,z,t)
The resultant deformation, (1,, is expressed as,
d,(x,y,z,t)="d(x,y,z,t)” (5.35)

Now that the program has calculated total piston deformation, the FEMstrain
subroutine is called. Here the strains for each element are calculated, as well as the
stresses at each node. The strain vector is calculated according to Eq. (2.93), and the
stress vector according to Eq. (2.96). Some further calculations are performed within the

F EMstrain subroutine. The equivalent or effective strain, seq, is calculated as,

seq = 2‘ ’51—:38—2— (5.36)

where,
e, --;:[(£x—£a)2 +(8y—ea)2 +(e, —€a)2:|
82 =§[(rxy)2 +(r... )2 +(ryz )2] (5.37)
a, =§(.,+gy+g,)

The von Mises stress component, own, is calculated from the stress components as,

101

0,0,, =[32—[(0x -ay)2 +(0x ‘02)2 +(0'y ‘02)2]

(5.38)
+ 3[(0xy )2 + (0",: )2 + (0y: )2 Hm

The principal stresses, 0,, 02, and 03 are given by the eigenvalues of the stress tensor:

0' 1 0 0 0x ny 0x2
0 0'2 0 = eig 0",), 0y ayz (5.39)
0 0 (73 0'12 6y: Oz

The stress intensity, mm, is calculated from the difference of the maximum and
minimum principal stresses. Given that a, > 02 > 03 then,

Jim 2 0'1 — 03 (5.40)

After all these calculations are made, the data is saved and the outer loop moves

to the next crank angle. The final converged solution of the hydrodynamic and contact

pressures will be periodic and will not depend on the initial guess for the total pressure.

47:
Ph(t)—Ph [t+;),
(5.41)

P

.(t)=a(r+i§)

It is observed that a periodic solution is usually achieved after the second cycle.
Convergence of pressure, however, is greatly dependent on clearance. A high clearance

will tend to converge up to ten times faster than a low clearance.

102

5.4 1r-fea Verification

As n—fea is a newly developed program, its predictions for both, the temperature
and the deformation are compared with COSMOSDesignSTAR’s predictions for
verification. A piston is considered with the thermal loads of in Table 4.2. A 100 kPa
uniform pressure is applied to the crown and the skirt, Also a body load of 9.81 ms'2 is
applied.

Figure 5.5 shows the nodal temperature as predicted by n-fea and
COSMOSDesignSTAR. The two predictions are in excellent agreement.

Figure 5.6 shows the prediction for the nodal resultant deformation. This
deformation is caused by both the thermal and mechanical loads. For the u-fea prediction
the contribution to the deformation due to the skirt loading is calculated via the Guyan

reduction technique (Section 2.4.5). The two predictions are also in excellent agreement.

 

 

 

 

 

330, .~ 7 7%
n-fea ‘
g 1,—:CO$A/LOSD-‘ . . 1 1
93 325» ' l
D ‘ 1 l
E ‘ 1
g 1 l l .
d) 1 1
'— 1
315‘ ,7,#,,, . i . ,. .-
1 201 401 601 801 1001 1201 1401
_ Node no.
E 325
(D
5
g 320
a)
‘e‘
31 n l l 1 l
11’ 340 850 860 870 880 890 900

Node no.

103

Figure 5.5: Nodal temperature as predicted by n— ea
and COSMOSDesignSTAR

40fi i 71 77 777‘ if 7 il;7::r::iflirri fiI
\ n-fea ‘ ‘
1; 43,98 MQSJZ-J .

 

OD
. 0

Res. Def. [p.m]
N
4 o!

 

 

 

 

ml

0‘ ,

1 20174017 601 801' 1001 714201 12101

Nodeno.

E40
EZOW
05
a) 0 J l 1 r r
‘1 340 850 860 870 880 890 900

Node no.

Figure 5.6: Nodal resultant deformation as predicted by 1t- ea
and COSMOSDesignSTAR

104

CHAPTER 6

RESULTS

6.1 Introduction

This chapter presents results that can be obtained from n- ea. It should be recalled
that this version of n—fea only performs a finite element analysis of the piston to calculate
piston deformation and piston skirt loading, assuming _r_1_o_ secondary motion of the piston.
Consequently the results are not typical for a real operating piston as secondary dynamics

greatly influence piston behavior.

6.2 Engine geometry and operating conditions
The engine geometry used for this simulation is the optical engine of Chapter 4.
Table 6.1 shows the required inputs for the engine geometry and operating conditions.

Figure 6.1 shows the pressure trace used for the simulation.

Table 6.1: Engine geometry and operating conditions

 

 

Parameter Value
Engine speed 1000 RPM
Bore diameter Simulation I: 90.2 mm
Simulation II: 90.1 mm
Stroke 90.6 mm
Connecting rod length 133 mm
Oil viscosity 0.02 Pa-s
Friction coefficient 0.15

 

 

105

 

30 | f T T f 7 l
' : : : l

Pressure [bar]

 

 

 

 

 

O 90 180 270 360 450 540 630 720
Crank angle [deg.]

Figure 6.1: Pressure trace used for simulation

6.3 Piston geometry and properties
The piston used for this simulation is a generic piston, designed to comply with
the engine geometry of Section 6.2. Table 6.2 shows the piston properties and Figure 6.2

shows the piston CAD model used for the simulation.

Table 6.2: Piston properties

 

 

Parameter Value
Mass 0.4536 kg
Diameter 90 mm
Pin offset —0.4 mm
Distance fi'om top of skirt to centre of pin 14.5 mm
Skirt waviness 3.50 um

 

 

106

 

Figure 6.2: Piston CAD model, (a) isometric view, anti-thrust side shown,
(b) underside, (c) front view, anti-thrust side shown,
(d) side view, red: thrust side

The material used for the piston is aluminum alloy 1345. Its material properties

are shown in Table 6.3.

Table 6.3: Material properties for aluminum alloy 1345

 

 

Propgrty Value
Modulus of Elasticity 69 GPa
Tensile yield strength 82.72 MPa
Poisson’s ratio 0.33
Density 2700 1(ng
Coefficient of thermal expansion 2.4 x 10'5 K'1
Thermal conductivity 200 Wm'lK'l

 

 

107

The thermal boundary conditions for the piston are shown in Table 6.4. The

values are as described in Section 4.2.1.

Table 6.4: Thermal boundary conditions

 

 

 

 

Heat transfer Ambient
Piston area coefficient temperature
(Wm"K“) 6)
Crown 200 470
Ring-pack area 1000 314
Skirt area (thrust and anti-thrust sides) 700 310
Underside 250 303

 

6.4 Mesh Refinement

The area of interest on the piston is the skirt area. In order to choose an adequate
element size for the mesh a mesh refinement study was performed. Figure 6.3 shows the
imported piston meshed geometry into n- ea. The global element size is 9.2502 mm. The
skirt area though is meshed using an element size of 3 mm (Figure 6.4). Five different
meshes were tested as of Table 6.5, varying the element size on the skirt.

To assess the quality of each mesh the hydrodynamic force on thrust and anti-
thrust sides as well as the total friction force were checked. On the thrust (Figure 6.5) side
meshes I, II, and III are in good agreement. Meshes IV and V give a slightly higher force
at each crank angle, with the maximum difference being at about 3 N. This implies an
error of about 6%. On the anti-thrust side (Figure 6.6) only meshes I and II are in very
g00d agreement. The maximum difference in this case is about 6 N. The total
hydrodynamic friction is in very good agreement for all five cases. This implies that the

velocity of the piston, Eq. (5.18), dominates in the evaluation of friction. Considering the

108

above mentioned results and mesh sizes (Table 6.5), mesh III was chosen for this study as

it is only a demonstrative study and high accuracy is not required.

 

Figure 6.3: Imported piston meshed geometry, global element size: 9.2502 mm, skirt
element size: 3 mm

 

 

 

 

o o
xs [mm] xs [mm]
(a) (b)
Figure 6.4: Skirt mesh, (a) thrust side, (b) anti-thrust side

109

Table 6.5: Mesh sizes

 

 

 

 

 

 

 

 

 

Global element Skirt element No. of No. of
Mesh . .
srze (mm) srze (mm) elements nodes
1 9.2502 5 3808 1273
II 9.2502 4 4133 1398
III 9.2502 3 4997 1728
IV 9.2502 2.5 6097 2119
V 9.2502 2 8034 2791
60 , . T
. 3 5 : — '
. -— ll .-
—— Ill
-— IV
_ V .4

F... [N]

 

 

 

 

__________________________________________

 

 

 

 

 

0 90 180 270 360 450 540 630 720

Crank angle [deg.]

Figure 6.5: Hydrodynamic force on thrust side for meshes of Table 6.5

110

Fh'2 [N]

 

 

 

 

 

 

 

 

 

Crank angle [deg.]

Figure 6.6: Hydrodynamic force on anti—thrust side for meshes of Table 6.5

Fm [N1

 

  

 

........................

,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

 

 

eeeee ‘ +++++ . +++++ 1

 

0 90 100 m 300 450 540 630 720

Crank angle [deg.]

Figure 6.7 : Total hydrodynamic friction force for meshes of Table 6.5

111

6.5 Periodicity and convergence

As an initial guess for the hydrodynamic pressure is required to start the
simulation, it is evident that the solution afier one cycle will not represent the actual one.
Consequently the simulation needs to be run for more than one cycle until the solution is
periodic. It has been observed for all the cases tested, that periodicity occurs afier the
second cycle. This can be seen in Figure 6.8 and Figure 6.9 where the total hydrodynamic
force on the piston is shown for the first three cycles for a piston-to-bore clearance of 100

um and 50 um respectively.

 

 

 

 

 

 

 

 

 

 

 

 

 

0 90 180 270 360 450 540 600 720
Crank angle [deg.]

Figure 6.8: Hydrodynamic force over three cycles, C = 100 pm

112

 

Cycle 1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0
n
0 m
l.‘ i.‘ n u.
. _ _ m w
n n 0 __ 6
........... 5.5% C
A s, .w
. k 5
_ 0 c ,
,0
. a ..
w .3... - iv... a
0M .W S d 7:3" + . ..
:40 e r m u “+8.11% *0
..m. 6 .0“ m 1*? 1t. «5% LC- 6
e v a s 3
mm. o a e _
n e M R
3 a m m
k 1 11111111 ..
0 n f 2
7 m .c
2C m H 0
-1 e351: 8
0 m4 xv;#lfifi, 1
8 d rw‘o; 4’;
1 0 .. .cfifirf
0 H 1.. my... 1.... a ,
9 0 .. , ea 1
0.” 1 I it}. t
o X Ari 1‘1”? i0
0 m 0
Ma 22.2920 02 332.com
F

 

100 um

C

cycle,

113

Crank angle [deg.]

Figure 6.10: Hydrodynamic pressure convergence at 3rd

Iterations

 

 

No. of Iterations

 

 

 

.__--._-_-.___-_I._ ________________________________

Residual
0
U1

 

 

 

 

0 90 180:;2‘70’13400 450 540 63077720
Crank angle [deg.]

Figure 6.11: Hydrodynamic pressure convergence at 3rd cycle, C = 50 um

Figure 6.10 and Figure 6.11 show the convergence of pressure. The number of
iterations required and the residual are shown for clearance of 100 um and 50 um. For a
tolerance of 103, the hydrodynamic pressure for the case of 50 um converges about ten

times slower than the 100 um clearance case.

6.6 Simulation Results

In the interest of space not all simulation results can be shown here. Appendix B
shows all the outputs of n—fea for one crank angle for the simulation with clearance of
100 pm. The chosen crank angle is at 365 degrees where the combustion gas pressure is
still relatively high, and the motion of the piston allows for the buildup of hydrodynamic

lubrication pressure. Some selected results are discussed below.

114

 

6.6.1 Temperature

Figure 6.12 shows the temperature distribution of the piston. As discussed earlier,
given sufficient operating time of the engine, the piston is assumed to be at steady state
over the cycle as temperature fluctuations occur on the outer layer of the piston about 2

mm thick. As expected the temperature is higher on the crown because of the combustion

gases, and it decreases moving towards the bottom of the skirt. Most of the heat is lost at

 

the ring-pack area where the heat transfer coefficient is highest.

Temp

   

Figure 6.12: Temperature distribution

6.6.2 Deformation

Figure 6.13 shows the x-component of the piston deformation at 405 degrees for
both 100 um and 50 pm of clearance. It can be seen that the latter clearance causes a
higher pressure buildup, which pushes the piston skirt more towards the piston center axis

(Figure 6.13b) than the 100 um clearance, (Figure 6.13a).

115

   

(b)
Figure 6.13: Deformation at 405 deg: x-component, (a) C = 100 um, (b) C = 50 um

6.6.3 Strain

Figure 6.14 shows the equivalent strain at 405 degrees for the clearance of 100
um. The whole piston is in the low range up to 300 microns. The high strain occurs
around the pinhole. This is because the restraint for the static finite element analysis is
put at the pinhole. The nodes there cannot move relative to their neighboring ones. Thus a

high strain arises at that area.

 

Figure 6.14: Equivalent strain at 405 deg., C = 100 um

116

6.6.4 Stress

Since the stress is a function of strain, the area with high stress is around the
pinhole (Figure 6.15). Considering the maximum distortion energy theory and applying
the 3D von Mises yield criterion [53] to the piston (Figure 6.16), it can be seen that a few
nodes fall outside the cylinder. (The radius of the cylinder is defined by the yield strength
of the material.) According to the von Mises criterion failure would occur at the pinhole,

however in real operating conditions the pinhole surface slides on the piston pin.

von Mises
(Past 107

 

 

Figure 6.15: von Mises stress at 405 deg., C = 100 pm

117

 

8
x10

 

x 10° 02 (pa) '1 01 (Pa)

Figure 6.16: 3D von Mises yield criterion at 405 deg., C = 100 um

6.6.5 Lubrication Pressure

Lubrication or hydrodynamic pressure provides support for the piston within the
cylinder. Very high pressures are undesirable as they can deform the flexible piston skirt,
and also sharp pressure gradients will increase friction. In this section the lubrication
pressure as calculated by Eq. (2.114) is shown at two different crank angles, 315 degrees,
mid-stroke of compression, and 405 degrees, mid—stroke of expansion. At 315 degrees
(Figure 6.17 and Figure 6.18) the pressure is zero everywhere on the skirt. In the absence
of piston secondary motion the squeeze film is negligible, thus the physical wedge term
of the Reynolds equation dominates. However the velocity of the piston here is negative
causing a positive relative velocity, Eq. (2.115), and thus a negative pressure, where with
the application of the half-Sommerfeld condition it goes to zero. The opposite effect is

seen at 405 degrees where the velocity of the piston is positive.

118

-

 

Comparing Figure 6.19 and Figure 6.20 with Figure 6.21 and Figure 6.22
respectively, it can be seen that halving the clearance creates a pressure about one order

of magnitude higher.
A study of these pressure profiles over a cycle can help improve piston design and

thus piston performance.

‘1

0.5.

Pressure [MPa]

 

40

    
 

50
20

 
 

0

 

 

ys [mm] 0 '50 xs [mm]

Figure 6.17: Hydrodynamic pressure on thrust side at 315 deg., C = 100 um

0.5 a

Pressure [MPa]

 

40\

    
 

50

 

 

 

20

   

0

 

 

Vs [mm] 0 ’50 xs [mm]
Figure 6.18: Hydrodynamic pressure on anti-thrust side at 315 deg, C = 100 um

119

 

 

 

 

0

y, [mm] ° '50 xs [mm]

Figure 6.19: Hydrodynamic pressure on thrust side at 405 deg., C = 100 um

50.
40.

E30.

20~

a]

Pressure

109

 

OJ
40

    
 

50
20 0
VS [mm] 0 '50 xs [mm]

Figure 6.20: Hydrodynamic pressure on anti-thrust side at 405 deg., C = 100 pm

120

 

H400~
u

g- 300.

.—.400~

g- 300.

Pressure

Pressure

500 ~

200~
100~

 

0s
40

    
 
    
 

50

20 0

VS [mm] 0 '50 xs [mm]

Figure 6.21: Hydrodynamic pressure on thrust side at 405 deg., C = 50 um

500 .

200~
100~

 

 

50

20 0

Y3 [mm] 0 '50 xs [mm]

Figure 6.22: Hydrodynamic pressure on thrust side at 405 deg., C = 50 pm

121

 

6.6.6 Oil Film Thickness

Figure 6.23 and Figure 6.24 show the oil film thickness at 405 degrees for the
clearance of 100 um for the thrust and anti-thrust sides respectively. The simulation
assumes fully flooded lubrication. From these figures it can be deduced that the flexibility
of the skirt increases moving towards the bottom of it. Also it can be seen that the piston
has suffered thermal expansion as the oil film thickness is in the range of 75 um to 90

um. Similar trends are observed at all crank angles.

 

 

90 ~
'5'
.3 854
3
O
E 80 ~
15

75 a

40

50
20 0
VS [mm] 0 '50 xs [mm]

Figure 6.23: Oil film thickness on thrust side at 405 deg, C = 100 um

122

 

 

 

E‘
385
a
0
C
i 80~
IE
75.
40
so
20 o
0 -50
mm
vs! 1 xslmm]

Figure 6.24: Oil film thickness on anti-thrust side at 405 deg, C = 100 um

6.6.7 Forces and Moments

Figure 6.25 shows the total hydrodynamic force experienced by the piston over a
cycle for both clearances, 100 um and 50 pm. As expected, the low clearance creates
much higher force. The hydrodynamic force, Eq. (5.16), is negative for both cases,
implying that the hydrodynamic force on the anti-thrust is higher. This is caused from the
thermal strains in conjunction with the pin offset. As described earlier, the pin offset is
towards the thrust side. Splitting the piston by a plane perpendicular to the pinhole axis
and coincident with the piston center axis (Figure 6.26), the anti-thrust side half has
slightly higher material volume. This allows for a higher thermal expansion on the anti—
thrust side half, as the pinhole is restrained and thermal expansion occurs radially away

fi'om the pinhole.

123

 

 

 

180 270 360 450 540 630 720

    

90

 

 

 

Crank angle [deg.]

Total hydrodynamic force on piston

Figure 6.25

 

 

Centre Axis Plane

 

Plane parallel to pinhole axis and coincident with piston axis

Figure 6.26

124

Thrust side Anti-thrust side

 

 

  

 

 

 

 

 

 

 

 

 

 

 

 

60 .
6
c 40
(D
3
9 20h
LI.
0 0 '
12 14 16 18 20 22 24 12 14 16 18 20 22 24
Res. deform. [um] Res. deform. [1110]
60 1 , 1 . .

i 1 ' Thrust side
Q40--------5.-------_5 ........ L-Anti-thrustside_
93 20 '

u.
0
12 14 16 18 20 22 24

Res. deform. [111“]

Figure 6.27: Histogram of resultant deformation on thrust and anti-thrust side due to
thermal expansion only

Figure 6.27 considers the histograms of the resultant deformation on thrust and
anti- thrust sides due to thermal expansion. The deformation of the anti—thrust side is
shifted to the right compared to the deformation of the thrust side. Also the anti-thrust
side deformation has higher frequency in the high values. This slightly higher
deformation on the anti-thrust side leads to a slightly smaller gap and consequently a
higher pressure buildup, which leads to a higher hydrodynamic force on the anti-thrust
side.

Also it can be seen from Figure 6.25 that the high clearance creates a force only
during the intake and expansion strokes. This is not the case for the low clearance. It
creates a force for a few degrees after the beginning of the compression stroke as well as

of the exhaust stroke. This is happens because of the squeeze film effect (second term on

125

 

the right hand side of Reynolds equation). The high pressures have created enough
deformation on the flexible skirt during the down-strokes, that at the beginning of the two
up-strokes where the physical wedge has no effect, the squeeze film builds up pressure,
as the piston skirt relaxes.

Figure 6.28 shows the total contact force. For these two cases no contact force
occurs. This is expected in the absence of piston secondary motion, unless the piston

experiences extremely high temperatures that will cause very high thermal expansion.

 

 

 

 

 

.—-C=100pm

05 ----- :———-i—---: ----- :———-[-—C=50um
Z... 0 ' I i I I l i
11.0 I i i i i :

 

 

l l

0 90 180 270 360 450 540 630 720
Crank angle [deg.]

 

Figure 6.28: Total contact force on piston

 

 

 

 

 

 

 

 

l l

 

Crank angle [deg.]

Figure 6.29: Total friction force on piston

126

 

 

Figure 6.29 shows the total friction force on the piston, Eq. (5.32). The low
clearance causes a higher friction, as expected. It is evident from the figure that the piston
velocity term dominates the value of friction. The friction profile looks like an inverted
piston velocity profile (Figure 3.5), but of different magnitude. Lubrication pressure
gradient has little effect on friction.

Figure 6.30 shows the total moment, Eq. (5.33), experienced by the piston. The
moment is relatively small for both cases, as it depends on the piston skirt centroid
relative to the pinhole center. Comparing the two, the low clearance creates a much

higher moment on the piston.

 

0.04
0.03
0.02
0.01

 

 

 

[Nm}

2 -0.01
-0.02 ~
-0.03 -

 

-0.04

 

 

 

 

 

 

 

90 180 2%0 360 450 540 630 720
Crank angle [deg.]

-0.05 ‘
0

Figure 6.30: Total moment on piston

127

 

 

CHAPTER 7

CONCLUSIONS AND RECOMMENDATIONS

7.1 Conclusions

In this thesis two studies were conducted: an evaluation of the parameterized
piston modeling method and development of a finite element program for the hill piston
model analysis.

The first study has demonstrated the advantage of using a parameterized piston
model for evaluating the piston characteristics, elastohydrodynamic lubrication,
dynamics, and friction. A careful parameterization of the full piston model lead to a
reduced model which provided reasonable simulation results, with a significant reduction
in computation time.

Parameterized piston analysis can be used as a guideline for optimal piston
design. Optimization can be done on the parameterized piston automatically where design
characteristics such as the skirt profile, the pin location, and the piston ovality can be
varied to reconstruct a new model without the need of external solid modeling software.
Computational time constraints would prohibit automatic optimization using a very
detailed FE model.

Another advantage of the parameterized piston method is that it requires minimal
FEM knowledge from the end user. The mesh is generated automatically. Also the
convective boundary conditions, the pressure loads, the body load, and the restraints are

applied automatically. Setting up a full model for a FE analysis can be a tedious job as

128

 

boundary conditions have to be applied to faces. For the piston presented in Chapter 4 the
fill] model has 782 faces whereas the parameterized model only 29.

The piston skirt has a strong influence on piston performance, such as lubrication,
dynamics and friction. It has been shown here that the parameterized piston demonstrates
a qualitative good agreement with the full model at the skirt area. Under thermal loading
they both deform away from the piston axis. Under mechanical loading they both deform
towards the piston center axis. Under the combined effect of mechanical and thermal
loading, they both deform away fiom the piston center axis.

There is a very good agreement between the effects of mechanical loading on the
two models, both qualitatively and quantitatively. The qualitative agreement is
maintained under the non-uniform skirt loading conditions. When the thermal expansion
is combined with the mechanical loading, the deformation results exhibit more of
qualitative agreement than quantitative. This disagreement on nodal deformation is
caused by the thermal expansion. As the parameterized model has a smaller volume,
125505 mm3, compared to the full model, 130833 mm3, its thermal expansion is less than
the full piston model. In the thermal analysis the piston is assumed to be at steady state
over the entire cycle, the in—cylinder air temperature is averaged, and mean values are
chosen for the heat transfer coefficients. Considering these assumptions the discrepancy
caused by the thermal expansion can be considered acceptable. The goal to obtain the
general behavior of the piston under both thermal loading and mechanical loading was
achieved here.

In the CASE analysis the deformation due to in-cylinder gas pressure, lubrication

pressure, and body load is evaluated at every crank angle. In this case the parameterized

129

 

 

piston demonstrates the same behavior as the fill] model. Thus the simulated
characteristics, like skirt profile and piston pin offset, are considered reasonable and can
be applied to the full model.

The second study has concentrated on the development of a finite element
program, specially tailored to conduct a piston analysis over a full four-stroke cycle. The
program couples heat transfer, solid mechanics, fluid mechanics and tribology under the
umbrella of finite elements. The program is able to import the meshed geometry of a
piston and perform a full cyclic piston analysis to predict the piston’s temperature
distribution, deformation, strains, and stresses. It also predicts skirt lubrication pressures,
oil film thicknesses, and forces and moments.

All these results are not typical of a piston at real operating conditions. A major
assumption was made, that the piston is moving along the center of the cylinder bore with
no freedom to move transversely or tilt. Also lubricant viscosity was assumed to be
constant along the length of the cylinder bore, and the cylinder bore was assumed rigid.
Despite these assumptions some insights were obtained for the piston modeling approach
and its behavior over a cycle.

The mesh size plays an important role in the accuracy of predicted results.
Although linear tetrahedral elements were used, their size can be controlled to minimize
error. For the assessment of the flexible skirt characteristics (forces and moments), a
mesh size of 2.5 mm would be adequate as the difference from the 2 mm element size is
very small. In the interest of computation time and since this was a demonstrative

simulation an element size of 3 mm was chosen for the skirt area. The rest of the piston

130

 

 

can be meshed with a bigger element size as no crucial information is extracted fiom it
consequently saving on computation time.

The internal combustion engine is a cyclic device. Thus a cyclic behavior is
expected for the predicted parameters. As piston modeling is treated as an initial value
problem, guessing an initial lubrication pressure, the simulation needs to be run for more
than one cycle to fade out the error caused by the initial guess. It has been observed that
cyclic behavior occurs afier the second cycle. Thus a simulation needs to be run for at
least two cycles for the predicted results to be considered valid.

Convergence of hydrodynamic pressure is directly affected by the nominal piston
to cylinder bore clearance. As the clearance decreases, the pressure requires more
iterations to converge. For all the cases studied the pressure never failed to converge. For
some crank angles though, this required up to 200 iterations.

A small clearance creates a higher pressure, as expected. In such a case the piston
skirt experiences less deformation. A small clearance, however, creates higher friction,
which is undesirable. At this stage, in the absence of secondary piston motion, no
conclusion can be made for the clearance size. From other efforts it has been shown that a
high clearance causes piston slap, and a too small clearance high friction. Consequently
there is a tradeoff between the two, as the hydrodynamic pressure is the one that supports
the piston within the cylinder bore.

The simulations predict that the area around the pinhole is the most highly
stressed area of the piston. The von Mises yield criterion even predicts failure at that area.

However, in real operating conditions, this prediction does not hold as the pinhole area is

131

 

 

not restrained but slides on the piston pin. For the simulation it has to be restrained to
provide sufficient boundary conditions for the solution of the static problem.

Finally it has been predicted that the piston pin offset affects lubrication pressure
and thus hydrodynamic force. Due to the offset the skirt face furthest from the pinhole, in
this case the anti-thrust side, experiences higher thermal expansion, creating a smaller

gap and thus building up a higher pressure.

7.2 Recommendations

The parameterized piston model achieves its goal. It predicts the piston’s general
behavior under the studied loads. The parameterized model can be improved by
incorporating the ring grooves, consequently bringing it closer to the full piston model.

1r-fea can currently predict a majority of the parameters affecting piston
performance. However, it needs further development before it can be confidently used in
a real simulation. The present developments form the backbone of a reliable piston
modeling program.

Secondary piston motion should be included in the program. The equations of
secondary motion of the piston and the equations of motion of the connecting rod should
be coupled with the deformation results and the Reynolds equation to predict piston
secondary motion within the cylinder bore.

Also, an algorithm should be developed to monitor the position of the piston land
areas, the areas between the ring grooves,. In the case of contact of any land with the

cylinder bore, the secondary piston dynamics would be altered.

132

 

The Reynolds equation also needs improvement. The modified Reynolds equation
developed by Patir and Cheng (1978) should be used to account for the effects of surface
texture on hydrodynamic lubrication.

The half-Sommerfeld boundary condition for the hydrodynamic pressure used in
this model violates mass conservation. The cavitation algorithm [17] should be used
instead. It solves for the Reynolds equation for density variations and then calculates
pressure, in the meantime accounting for mass conservation.

The asperity contact model should be improved as well. Currently the formula
developed by Zhu et al. [59] for the contact of an aluminum piston with a cast iron
cylinder bore is used. This limits the model to just these two materials. The Greenwood-
Tripp model [19], which accounts for the contact of two nominally flat rough surfaces, is
recommended.

A flexible cylinder bore should also be included in the model, as in real operating
conditions the combustion gas pressure, lubrication pressure, and temperature gradients
deform the cylinder bore. Finally the lubricant viscosity should be modeled to vary with
temperature as this is the case with conventional lubricants. The cylinder bore is at
different temperatures along its length as well as the piston skirt. These affect the
lubricant viscosity in real operating conditions.

In conclusion, the parameterized piston model demonstrates a similar thermal and
mechanical behavior as the full piston model. The 1r-fea program can import the meshed
geometry of a full piston model and perform an analysis over a four-stroke cycle to
predict piston deformation, strains and stresses, as well as hydrodynamic and contact

forces.

133

 

 

APPENDIX A

134

APPENDIX A

FLOWCHARTS OF THE n—fea PROGRAM AND OF ITS SUBROUTINES

Start n—fea

Load the mesh data files

Load the piston material properties,
the engine geometry,
the boundary conditions, and the pressure trace

Calculate the piston motion

GO to FEMmeshread

GO to FEMthermaI

GO to FEMstatic

Sort nodes on skirt

GO to FEMGuyan

Map the skirt profile to a 2-D

coordinate system and renumber
the skirt connectivity matrix

 

Figure A.l: Flow chart for the main program of n-fea (Part 1 of 2)

135

V
I
I
I

    

For all the crank angle degrees

Calculate the deformation due to gas pressure

 

Calculate the deformation due to piston acceleration

WHILE the pressure residual greater than 10 '3 .

-'
I
l
I

GO to FEMs‘taticLUB

Calculate the total deformation, gap,
squeeze film, and contact pressure

GO to ReEq

Calculate the total pressure on skirt

Calculate the hydrodynamic, contact,
and fiiction forces

 

Calculate the piston deformation due to the skirt loading
and the total piston deformation

GO to FEMstrain

Save data

III-I-I-I-I-

End u-fea

Figure A.2: Flow chart for the main program of n—fea (Part 2 of 2)

136

Enter from 1r- ea

Load the thermal mesh data file

FOR all the lines

 

Read the nodal coordinates

Read the connectivity matrix
and the faces to be plotted

Read the convective boundary
conditions

Load the static mesh data file

FOR all the lines

Read the restraint nodes

Read the uniform pressure
boundary conditions

Read the body load

 

Exit to n—fea

Figure A.3: Flowchart for FEMmeshread subroutine

137

Enter from 1r-fea

Initialize arrays and matrices

FOR all the elements

Calculate the element volume and characteristics

 

Calculate the element stiffness matrix

 

. - as u - - -- FOR all the faces

IF convective boundary

Calculate the element load vector and matrix

 

Assemble the global stiffness matrix and load vector

Solve for the temperature
using conjugate gradient method

Exit to n—fea

 

Figure A.4: Flowchart for FEMthermal subroutine

138

»------------------------------------.

- Enter from n-fea

Initialize arrays and matrices

FOR all the elements

Load the element characteristics

Calculate the element stiffness matrix

 

Calculate the body force load vector

Calculate the thermal strain load vector

FOR all the faces

 

.‘ \

I.
I
I
I
I

IF uniform pressure boundary

—’

IF crown flag

Calculate the surface force load vector

   

...,f _“- IF thrust side flag

 

  
  

I'

I

1

: Build the thrust side skirt connectivity matrix

I

.” m. ‘ 1F anti-thrust side flag

1

I

I

: Build the anti-thrust side skirt connectivity matrix

it

Figure A.5: Flowchart for FEMstatic subroutine (Part 1 of 2)

139

“a

Assemble the global stiffness

Assemble the global thermal stress load vector

Assemble the global body force load vector

Assemble the global surface force load vector

-_-_-___..--..__--———-—_-—-___-i

Solve for the deformation due to thermal strains
using conjugate gradient method

Solve for the deformation due to body force
using conjugate gradient method

Solve for the deformation due to surface forces
using conjugate gradient method

Exit to n—fea

 

Figure A.6: Flowchart for FEMstatic subroutine (Part 2 of 2)

140

 

 

 

Enter from n-fea

FOR all the nodes on the skirt

Calculate the reduced stiffness matrix

Calculate the coordinate transformation matrix

Calculate the skirt compliance matrix

Exit to 1r- fea

Figure A.7: Flowchart for FEMGuyan subroutine

141

 

Enter from n- ea

Initialize arrays and matrices

FOR all the elements that have faces on the skirt

Load the element volume and characteristics

FOR all the faces

IF uniform pressure boundary condition

IF thrust OR anti-thrust flag

Calculate the load vector due to hydrodynamic pressure

Solve for the skirt deformation

Exit to Ir-fea

 

Figure A.8: Flowchart for FEMstaticL UB subroutine

142

 

 

'-—-u----————-—---—----—-—

 

 

Enter from n- ea

Initialize arrays and matrices

FOR all the elements

Load the element area and characteristics

 

Calculate the element stiffiiess matrix

Calculate the element load vector
due to the physical wedge

Calculate the element load vector
due to squeeze film

Assemble the global stiffness matrix
and global load vector

IF global stiffness matrix IS sparse

Solve for the hydrodynamic pressure
using conjugate gradient method

 

IF global stiffness matrix IS NOT sparse

Solve for hydrodynamic pressure
by inversion

Exit to n—fea

Figure A.9: Flowchart for ReEq subroutine

143

Enter from u—fea

Initialize arrays and matrices

FOR all the elements

Load the element volume and characteristics

Calculate the strain vector from known deformation

Calculate the stress vector from strain vector

Calculate the equivalent strain

Calculate the von Mises and principal stresses
and stress intensity

Exit to u-fea

 

Figure A.10: Flowchart for FEMstrain subroutine

144

 

 

 

 

APPENDIX B

 

145

 

APPENDIX B

SIMULATION RESULTS

The simulation results shown in this appendix are for the engine and piston
described in Chapter 6. All the results shown are at 365 crank angle degrees, 5 degrees
alter the expansion stroke. This crank angle was chosen because the combustion gas
pressure is still very high and the motion of the piston allows for hydrodynamic pressure
to be built. All these results are available for each crank angle, which make it feasible to
create animations to visualize piston deformation, stress, lubrication pressure or 011 film

thickness over a cycle.

' z in.»

'1» ' . .r .
.‘:1.n*?~.,‘r«mou-.w——~tfi.. ,.

 

Figure 8.1: Temperature distribution

146

 

x-dlm 5
("0x 10'

 

   

Figure 8.2: Deformation: x-component

y-dlm 5
(m)x 10'
2

 

   

Figure 8.3: Deformation: y-component

147

 

 

 

Figure 3.4: Deformation: z-component

Res
("0x 10-5

   

Figure 3.5: Deformation: resultant

148

 

 

Figure 3.7: Normal strain: x-component

149

 

Figure B.8: Normal strain: y-component

 

Figure 3.9: Normal strain: z-component

150

 

Figure 8.11: Shear strain: xz-component

151

 

Figure 8.13: Equivalent strain

152

 

 

Figure 3.14: Normal stress: x-component

0'

y 7
(P331 10

   

Figure 3.15: Normal stress: y-component

153

 

 

 

 

Figure B.16: Normal stress: z-component

1:

xy 7
(Pa): 10

   

Figure B.17: Shear stress: any-component

154

 

 

Figure 8.18: Shear stress: xz-component

‘11

112 7
(Pair 10

 

Figure B.19: Shear stress: yz-component

155

 

 

 

Figure B.20: Principle stress: lSt principal direction

°2

(”it 107

 

Figure B.21: Principle stress: 2nd principal direction

156

 

 

Figure B.22: Principle stress: 3rd principal direction

von Mlses

(Pair 1o7

    

Figure 8.23: von Mises stress

157

 

 

Figure 8.24: Stress intensity

 

 

 

 

 

 

 

02 (Pa) 0'1 (Pa)

Figure B.25: 3D von Mises yield criterion

158

 

 

10~
8.
6‘

Pressure [kPa]

2x

 

0s
40

  

50

     
  

20 0

Y9 [mm] 0 :50 x8 [mm]

Figure 3.26: Hydrodynamic pressure: thrust side

10~
8t
6.

Pressure [kPa]

2.

 

0 a
40
50

 
   
 

20 0

Y5 [mm] 0 “'50 xs [mm]

Figure 3.27: Hydrodynamic pressure: anti-thrust side

159

Thickness [p.fl'l]
CD
0|

 

hfl
DO!
U

0
U'l

 

nu
(DUI
I!

n} 3,»... ”I ' 4 “‘51:;
Fhfieasm‘.‘ -

 

 

50

2O 0

Vs [mm] o '50 X. [mm]
Figure 3.28: Oil film thickness: thrust side
it“ “I “win is»...
Seenturtfiv“
50

20 0

Y5 [Mm] 0 '50 XS [mm]

Figure B.29: Oil film thickness: anti-thrust side

160

REFERENCES

161

 

10.

ll.

12.

13.

REFERENCES

. Annand, W. J. D., 1963, “Heat Transfer in the Cylinders of Reciprocating Internal

Combustion Engines,” Proceedings of the Institution of Mechanical Engineers, Vol.
177, pp. 973-990.

Beer, F. P., Johnston, R. E., DeWolf, J. T., Mechanics of Materials, McGraw-Hill
Higher Education, New York, 2002.

Bejan, A., Advanced Engineering Thermodynamics, John Wiley & Sons, Inc., 1997.
CASE Theoretical Manual, 2005, Mid Michigan Research, Okemos, Michigan.

Cheney, W., Kincaid, D., Numerical Mathematics and Computing, Brooks/Cole
Publishing Company, New York, 1999.

Chui, B.K., Schock, H.J., 2002, “Computational Analysis of Piston-Ring Wear and
Oil Consumption for An Internal Combustion Engine,” ASME ICE, Paper No.:
ICEF2002-529.

Chui, B.K., Schock, H.J., Fedewa, A.M., Richardson, DE, Shaw, T., 2005, “Coupled
Model of Partially Flooded Lubrication and Oil Vaporization in an Internal
Combustion Engine”, ASME ICE, Paper No.: ICE82005-1077.

COSMOSM 2.95 Electronic Documentation, Structural Research and Analysis
Corporation, 2005.

Dursunkaya, Z., Keribar, R., 1992, “Simulation of Secondary Dynamics of
Articulated and Conventional Piston Assemblies,” SAE Paper 920484.

Dursunkaya, Z., Keribar, R., Ganapathy, V., 1994, “A Model of Piston Secondary
Motion and Elastohydrodynamic Skirt Lubrication,” Journal of T ribology, Vol. 116,
pp. 777-785.

Duyar, M., Bell, D., Perchanok, M., 2005, “A Comprehensive Piston Skirt
Lubrication Model Using a Mass Conserving EHL Algorithm,” SAE Paper 2005-01-
1640.

Ejakov, M.A., 2001, “Modeling of Axial and Circumferential Ring Pack
Lubrication,” ASME ICE, Paper No.: 2001-ICE-433.

Ejakov, M.A., Diaz, AR, and Schock, H.J., 1999, “Numerical Optimization of Ring-
Pack Behavior,” SAE Paper 1999-01—1521.

162

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

26.

27.

28.

Ejakov, M.A., Schock, H.J and Brombolich, L.J., 1998, “Ring Twist Modeling for an
IC Engine,” SAE Paper 982693.

Ejakov, M.A., Schock, H.J., Brombolich, L.J., Carlstrom, CM. and Williams, R.L.,
1997, “Simulation Analysis of Inter-Ring Gas Pressure and Ring Dynamics and Their
Effect on Blow-by,” ASME ICE, Paper No.: 97-ICE-50.

Ejakov, M.A., Yeager, DA, and Kayupov, M., 2003, “Reliability and Robustness
Analysis of Engine Ring-pack Performance,” SAE Paper 2003-01-1221.

Elrod, H. G., 1981, “A Cavitation Algorithm,” Journal of T ribology, Vol. 103, pp.
350-354.

Goenka, P. K. and Meemik, P. R., 1992, “Lubrication Analysis of Piston Skirts,” SAE
Paper 920490.

Greenwood, J. A., Tripp, J. H., 1971, “The Contact of Two Nominally Flat Rough
Surfaces,” Proceedings of the Institution of Mechanical Engineers, Vol. 185, pp. 625-
633.

Guyan, R. J ., 1965, “Reduction of Stiffness and Mass Matrices,” AIAA Journal, Vol.
3, No. 2, p. 380.

Guzella, L., Onder, C. H., Introduction to Modeling and Control of Internal
Combustion Engine Systems, Springer, New York, 2004.

Harrison, H. R., Nettleton, T., Advanced Engineering Dynamics, John Wiley & Sons,
Inc., New York, 1997.

Heinrich, J. C., Pepper, D. W., Intermediate Finite Element Method: Fluid Flow and
Heat Transfer Applications, Taylor & Francis, Philadelphia, 1999.

Herbst, M. H., Priebsch, H. H., 2000, “Simulation of Piston Ring Dynamics and Their
Effect on Oil Consumption,” SAE Paper 2000-01-0919.

Heywood, J. B., Internal Combustion Engine Fundamentals, McGraw Hill, Inc., New
York, 1988.

Hibbeler, R. C., Engineering Mechanics: Dynamics, Prentice Hall PTR, Upper Saddle
River, 2001.

Incropera, F. P., DeWitt, D. P., Introduction to Heat Transfer, John Wiley & Sons,
Inc., New York, 2002.

Kaliakin, V. N., Introduction to Approximate Solution Techniques, Numerical
Modeling, and Finite Element Methods, Marcel Dekker, Inc., New York, 2002.

163

 

29.

30.

31.

32.

33.

34.

35.

36.

37.

38.

39.

40.

41.

42.

43.

Kays, W. M., Convective Heat and Mass Transfer, McGraw-Hill Higher Education,
Boston, 2005.

Keribar, R., Dursunkaya, Z., 1992, “A Comprehensive Model of Piston Skirt
Lubrication,” SAE Paper 920483.

Keribar, R., Dursunkaya, Z., Ganapathy, V., 1993, “An Integrated Design Analysis
Methodology to Address Piston Tribological Issues,” SAE Paper 930793.

Khonsari, M. M., Booser, E. R., Applied T ribology: Bearing Design and Lubrication,
John Wiley & Sons, Inc., New York, 2001.

Kiusalaas, J., Numerical Methods in Engineering with MATLAB, Cambridge
University Press, Cambridge, 2005.

Kreyszig, B, “Advanced Engineering Mathematics,” John Wiley & Sons, Inc., New
York, 2005.

Li, C. H., 1982, “Piston Thermal Deformation and Friction Considerations,” SAE
Paper 820086.

Li, C. H., 1986, “Thermoelastic Behavior of an Aluminum Diesel Engine Piston,”
SAE Paper 860136.

Li, D. F., Rohde, S. M., Ezzat, H. A., 1982, “An Automotive Piston Lubrication
Model,” ASLE Transactions, Vol. 26, pp. 151-160.

Ludema, K. C., Friction, Wear, Lubrication: A Textbook in T ribologz, CRC Press,
Inc., New York, 1996.

Moran, M. J., Shapiro, H. N., Fundamentals of Engineering Thermodynamics, John
Wiley & Sons, Inc., New York, 2004.

Oh, K. R, Li, C. H. and Goenka, P. K., 1987, “Elastohydrodynamic Lubrication of
Piston Skirts,” Journal of T ribology, Vol. 109, pp. 356-362.

Oppenheim, A. K., Combustion in Piston Engines: Its Technology, Evolution,
Diagnosis, and Control, New York, Springer, 2004.

Panayi, A., Schock, H., Chui, B.K., Ejakov, M., 2006, “Parameterization and FEA
Approach for the Assessment of Piston Characteristics,” SAE Paper 2006-01-429.

Panton, R. L., Incompressible Flow, John Wiley & Sons, Inc., New York, 1996.

164

 

44.

45.

46.

47.

48.

49.

50.

51.

52.

53.

54.

55.

56.

57.

58.

Pepper, D. W., Heinrich, J. C., The Finite Element Method: Basic Concepts and
Applications, Hemisphere Publication Corp., Washington, 1992.

Perchanok, M., 2000, “Modeling of Piston-Cylinder Lubrication with a Flexible Skirt
and Cylinder Wall,” SAE Paper 2000-01-2804.

Pope, S. B., Turbulent Flows, Cambridge University Press, New York, 2000.

Rao, S. 8., Applied Numerical Methods for Engineers and Scientists, Prentice Hall
PTR, Upper Saddle River, 2005.

Rao, S. S., The Finite Element Method in Engineering, Butterworth-Heinemann,
Boston, 1999.

Reddi, M. M., 1969, “Finite-Element Solution of the Incompressible Lubrication
Problem,” Journal of Lubrication Technology, Vol. 91, pp. 524-533.

Reddy, J. N., An Introduction to the Finite Element Method, McGraw-Hill, Inc., New
York, 1993.

Reddy, J. N., Gartling, D.K., The Finite Element Method in Heat Transfer and Fluid
Dynamics, CRC Press, Inc., Boca Raton, 1994.

Thompson, E. G., Introduction to the Finite Element Method, Theory, Programming,
and Applications, John Wiley & Sons, Inc., New York, 2005.

Ugural, A. C., Fenster, S. K., Advance Strength and Applied Elasticity, Prentice Hall
PTR, Englewood Cliffs, 2003.

Wong, V.W., Tian, T., Lang, H., Ryan, J.P., Sekiya, Y., Kobayashi, Y., Aoyama, S.,
1994, “A Numerical Model of Piston Secondary Motion and Piston Slap in Partially
Flooded Elastohydrodynamic Skirt Lubrication,” SAE Paper 940696.

Woschni, G., 1967, “Universally Applicable Equation for the Instantaneous Heat
Transfer Coefficient in the Internal Combustion Engine,” SAE Paper 670931.

Woschni, G., 1979, “Determination of Local Heat Transfer Coefficients at the Piston

of a High Speed Diesel Engine by Evaluation of Measured Temperature
Distribution,” SAE Paper 790834.

Wu, H., Chiu, G, 1986, “A Study of Temperature Distribution in a Diesel Piston —
Comparison of Analytical and Experimental Results,” SAE Paper 861278.

Zhu, D. Cheng, H.S., Arai, T. Hamai, K., 1992, “A Numerical Analysis for Piston

Skirts in Mixed Lubrication - Part I: Basic Modeling,” Journal of T ribology, Vol.
114, pp. 553-562.

165

 

59. Zhu, D. Hu, Y., Cheng, H.S., Arai, T., amai, K., 1993, “A Numerical Analysis for
Piston Skirts in Mixed Lubrication — Part II: Deformation Considerations,” Journal of
Tribology, Vol. 115, pp. 125-133.

166

 

 

I"Lilliiliifllj[[11111]]I