THREE-DIMENSIONAL MULTI-PHYSICS MODELING METHODOLOGY TO STUDY ENGINE CYLINDER-KIT ASSEMBLY TRIBOLOGY AND DESIGN CONSIDERATIONS By Sadiyah Sabah Chowdhury A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering - Doctor of Philosophy 2021 ABSTRACT THREE-DIMENSIONAL MULTI-PHYSICS MODELING METHODOLOGY TO STUDY ENGINE CYLINDER-KIT ASSEMBLY TRIBOLOGY AND DESIGN CONSIDERATIONS By Sadiyah Sabah Chowdhury Engine cylinder-kit tribology is pivotal to durability, emission management, friction, oil consumption, and efficiency of the internal combustion engine. The piston ring pack dynamics and the flow dynamics are critical to engine cylinder-kit tribology and design considerations. A three- dimensional (3D), multi-physics methodology is developed to investigate the liquid oil- combustion gas transport and oil evaporation mechanisms inside the whole domain of the cylinder kit assembly during the four-stroke cycle using multiple simulation tools and high-performance computing. First, a CASE (Cylinder-kit Analysis System for Engines) 1D model is developed to provide necessary boundary conditions for the subsequent steps of the chain of simulations. Next, the ring-bore and ring groove conformability along with the twist angle variation across the circumference are investigated by modeling a twisted ring via a 3D ring FEA contact model. The ring twist induces change in ring location which subsequently changes the cylinder kit geometry dynamically across the cycle. The dynamically varying geometries are generated using the LINCC (Linking CASE to CFD) program. Finally, a three-dimensional multiphase flow model is developed for the dynamic geometries across the cycle using CONVERGE. The methodology is first applied on a small-bore (50 mm) engine running at 2000 rpm. Next, a CASE 1-D model is developed and calibrated via HEEDS across a range of load-speed operating conditions of a Cummins 6-cylinder, 137.02 mm bore, Acadia engine. The 1800 RPM, full load condition with a positively twisted second ring is selected for the experimental validation of the 3-D methodology. A study of the second ring dynamics in the small-bore engine showed the effect of negative ring twist on the three-dimensional fluid flow physics. The oil (liquid oil and oil vapor) transport and combustion gas flow processes through the piston ring pack for the twisted and untwisted geometry configurations are compared. A comparison with the untwisted geometry for this cylinder-kit shows that the negatively twisted second ring resulted in a higher blowby but lower reverse blowby and oil consumption. The comparison of the model predicted oil consumption with existing literature shows that oil consumption is within the reasonable range for typical engines. The blowby, second land pressures and third land pressures comparison with the experimental results of Cummins Acadia engine showed considerable agreement. The reverse blowby and oil consumption along with the liquid oil and oil vapor mass fraction distribution pattern across the cycle are also analyzed. In the later section of this work surface texture characterization of a novel Abradable Powder Coating (APC) and stock piston skirt coatings of a Cummins 2.8 L Turbo engine is conducted. The surface texture and characteristic properties varying across the piston skirt are obtained and analyzed via a 3D optical profiler and OmniSurf3D software. The engine operating conditions are found through a combination of measurements, testing, and a calibrated GT-Power model. The variable surface properties along with other geometric, thermodynamic, material properties are utilized to build a model in CASE for both APC and stock coated pistons. The Surface texture analysis shows that the APC coating has a unique feature of mushroom cap-like surface and deeper valleys that could potentially be beneficial for lubrication and oil retention. Comparison of different performance parameters from CASE simulation results shows that APC has the potential to be a suitable candidate for piston skirt coating. ACKNOWLEDGEMENTS First, I would like to express my sincere gratitude to my advisor Dr. Harold Schock for being an amazingly inspirational mentor throughout my PhD journey. Thanks for pushing me to my greatest limits and helping me exceed bounds that I set for myself. I would also like to thank Dr. George Zhu, Dr. Farhad Zaberi, and Dr. Dirk Colbry for serving on my dissertation committee. Special thanks to Dr. Colbry for your guidance on high performance computing. I am thankful to Mid-Michigan Research LLC, Cummins Inc, Line2LineCoatings LLC and Commodo Materials LLC for supporting this work. Throughout my PhD years I have collaborated with cross functional teams of these companies. Being part of these projects and seeing the direct impact of the works in industry application has been a truly rewarding experience. My utmost gratitude to Dr. Larry Brombolich and Dr. Andreas Panayi for sharing their knowledge and experience in cylinder kit modeling. Scott Baeder, thanks for our hours of discussions on system configuration, software licensing, GUI implementation and what not. I am grateful to Dr. Chao Cheng and Dr. Ali Kharazmi- who started the three-dimensional modeling effort of the cylinder- kit. Their insights helped me understand the basic of the problems and the challenges to face. My sincerest gratitude to Ali for guiding me since the very beginning and being the key contact person for all the collaboration works with Cummins inc. I would like to thank Andy Suman for the opportunity to work on the novel abradable powder coatings. Our discussions on coating materials and modeling helped me dig deeper into the problem. I am thankful to Dr. Donald Cohen of Michigan Metrology and Dr. Mark Malburg of Digital Metrology for sharing their knowledge to help me understand the surface texture of the coatings. I would like to take this opportunity to acknowledge the support and friendship of my fellow EARL iv members. Shout out to the awesome Spartan undergrads who worked with me- Brian Fedewa, Helen Miller, Abhyuday Rastogi, Matthew Donahue and Daniel Nicklowitz. You guys are the best! My officemate Dr. Sedigheh (Ati) Tolou had been very supportive and welcoming since day one. Thanks, Ati! My deep appreciation goes to Tom Stuecken, Yidnekachew Ayele, Brian Rowley, John Przybyl and Brian Deimling for all the coffee breaks and making my time at EARL enjoyable. My greatest gratitude goes to my friends and family who endured this long process with me, always offering support and love. Thanks to Dr. Rebeka Sultana Lubna and Sadia Sayed for always having my back. Many thanks to My brother Panini Amin Chowdhury for believing in me. To Dr. Cyrus Ashok Arupratan Atis, thank you so much for your love, guidance, and immense support during tough times, and for rejoicing happy moments with me during good times. The two most important persons behind my story are my parents. My father Dr. Monzur-Ul-Amin Chowdhury instilled the idea of “never stop learning” since I have known life. My final words are for the person who would have been the happiest on earth seeing me achieve this- my mother Rokeya Begum. It has been five years since that moment when my life changed forever. Not a single day passed when I did not feel your presence giving me strength. Thanks for never letting me give up on anything I strive for, no matter how difficult the journey is. And, like my every little achievement, this is for you, Ammu! v TABLE OF CONTENTS LIST OF TABLES ......................................................................................................................... ix LIST OF FIGURES ....................................................................................................................... xi Chapter 1 INTRODUCTION .......................................................................................................... 1 1.1 Background and Motivation ............................................................................................... 1 1.2 Structure of Dissertation ..................................................................................................... 8 Chapter 2 THREE-DIMENSIONAL MULTI-PHASE PHYSICS-BASED MODELING METHODOLOGY TO STUDY ENGINE CYLINDER-KIT ASSEMBLY TRIBOLOGY AND DESIGN CONSIDERATIONS- PART I ....................................................................................... 9 2.1 Introduction ......................................................................................................................... 9 2.2 Workflow .......................................................................................................................... 10 2.2.1 One-Dimensional Simulation & Geometry Generation ........................................... 11 2.3 Three-Dimensional Multi-phase Model Development ..................................................... 13 2.3.1 Volume of Fluid (VOF) Model ................................................................................ 13 2.3.2 Governing Equations & Solver Setup ...................................................................... 16 2.3.3 Homogeneous Relaxation Model ............................................................................. 20 2.3.4 Boundary & Initial Conditions................................................................................. 22 2.4 Results and Discussion ..................................................................................................... 26 2.4.1 Grid Sensitivity Test ................................................................................................ 26 2.4.2 Circumferential Flow ............................................................................................... 29 2.4.3 Estimating Blow-by ................................................................................................. 31 2.4.4 Reverse Blow-by ...................................................................................................... 33 2.4.5 Pressure Distribution ................................................................................................ 35 2.4.6 Oil Transport ............................................................................................................ 35 2.4.7 Oil Evaporation ........................................................................................................ 46 2.5 Conclusions ....................................................................................................................... 49 Chapter 3 THE EFFECT OF RING-GROOVE GEOMETRY ON ENGINE CYLINDER-KIT ASSEMBLY USING THREE-DIMENSIONAL MULTIPHASE PHYSICS-BASED MODELING METHODOLOGY-PART II .................................................................................. 51 3.1 Abstract ............................................................................................................................. 51 3.2 Introduction ....................................................................................................................... 52 3.3 Workflow .......................................................................................................................... 55 3.4 Results & Discussion ........................................................................................................ 65 3.4.1 Circumferential flow ................................................................................................ 65 3.4.2 Estimating Blow-by ................................................................................................. 68 3.4.3 Estimating Reverse Blow-by & Oil Consumption from the Cylinder Kit ............... 70 3.4.4 Comparison of Oil Mass Fraction Distribution Along the Piston Length between twisted and untwisted geometry........................................................................................ 75 3.5 Conclusions ....................................................................................................................... 82 vi Chapter 4 ONE DIMENSIONAL OPTIMAZATION STUDY OF CUMMINS ACADIA ENGINE ....................................................................................................................................................... 84 4.1 Introduction ....................................................................................................................... 84 4.2 Methodology ..................................................................................................................... 86 4.2.1 Problem Setup in CASE........................................................................................... 86 4.2.2 Optimization Study Setup in HEEDS ...................................................................... 89 4.3 Results and Discussion ..................................................................................................... 91 4.3.1 Optimization at all operating conditions using ten variables. .................................. 91 4.3.2 Optimization at five operating conditions using the gas flow coefficient used for 1800 rpm 100% ................................................................................................................. 95 4.4 Summary and Conclusions ............................................................................................... 99 Chapter 5 EXPERIMENTAL VALIDATION OF THREE-DIMENSIONAL MULTIPHASE PHYSICS-BASED MODELING METHODOLOGY FOR ENGINE CYLINDER-KIT ASSEMBLY ............................................................................................................................... 101 5.1 Abstract ........................................................................................................................... 101 5.2 Introduction ..................................................................................................................... 102 5.3 Workflow ........................................................................................................................ 103 5.3.1 Block 1(One dimensional model in CASE) ........................................................... 103 5.3.2 Block 2 (Three-dimensional Finite Element Analysis (FEA) contact model) ....... 108 5.3.3 Block 3(Dynamically varying three-dimensional geometry generation)............... 116 5.3.4 Block 4 (Three-dimensional Computational Fluid Dynamics (CFD) model in CONVERGE ................................................................................................................... 117 5.4 Results & Discussion ...................................................................................................... 120 5.4.1 Circumferential flow .............................................................................................. 120 5.4.2 Blow-by.................................................................................................................. 122 5.4.3 Estimating Reverse Blow-by & Oil Consumption from the Cylinder Kit ............. 124 5.4.4 Land Pressure Comparison .................................................................................... 125 5.4.5 Liquid Oil and Oil Vapor Mass Fraction Distribution Along the Piston Length... 126 5.4.6 Speed Comparison: MSU HPC Vs. Argonne ........................................................ 130 5.5 Summary and Conclusions ............................................................................................. 134 Chapter 6 TRIBOLOGICAL PERFORMANCE ASSESSMENT OF ABRADABLE POWDER COATED PISTONS CONSIDERING PISTON SKIRT GEOMETRY AND SURFACE TOPOGRAPHY.......................................................................................................................... 135 6.1 Abstract ........................................................................................................................... 135 6.2 Introduction ..................................................................................................................... 136 6.3 Piston Surface Morphology Analysis ............................................................................. 140 6.4 Comparing Surface Properties of Worn Stock and Worn Coated (Post-run) ................. 147 6.4.1 Worn Coated vs Worn Stock at Side Location ...................................................... 148 6.4.2 Worn Coated vs Worn Stock at Center Location................................................... 149 6.5 CASE Model Setup ......................................................................................................... 152 6.5.1 Block 1(Input parameters and files generation) ..................................................... 152 6.5.2 Block 2 (Running model using CASE-PISTON) .................................................. 158 6.6 CASE Simulation Results ............................................................................................... 164 6.7 Summary and Conclusions ............................................................................................. 173 vii Chapter 7 Concluding Remarks and Recommendations for Future Work ................................. 176 7.1 Coated Piston Performance Assessment: Concluding Remarks and Recommendations for Future Work .................................................................................................................... 176 7.2 3D Multiphysics Methodology: Concluding Remarks and Recommendations for Future Work .......................................................................................................................... 178 APPENDIX ................................................................................................................................. 188 BIBLIOGRAPHY ....................................................................................................................... 193 viii LIST OF TABLES Table 2.1 Engine operating condition ........................................................................................... 12 Table 2.2 Convergence tolerance values ...................................................................................... 19 Table 2.3 Temperature boundary condition .................................................................................. 25 Table 2.4 Boundaries assigned to the regions ............................................................................... 25 Table 2.5 Simulation time and grid details of the cases ............................................................... 27 Table 3.1 Engine operating conditions [83] .................................................................................. 56 Table 3.2 Selected crank angles at which the ring location and OD change ................................ 62 Table 3.3 Selected specifications of the engines used in literature............................................... 73 Table 4.1 Range for variables ....................................................................................................... 89 Table 4.2 Best Design values of all ten variables for different operating conditions ................... 91 Table 4.3 Blowby and gas flow coefficients used across different operating conditions ............. 92 Table 4.4 Best design values of seven variables while using the same gas flow coefficients as in 1800 RPM 100% in the other operating conditions ..................................................... 96 Table 5.1 Engine operating conditions ....................................................................................... 104 Table 5.2 Main parameters for the ring ...................................................................................... 110 Table 5.3 Acadia Ring Tension .................................................................................................. 114 Table 5.4 Selected crank angles at which the ring location and OD change .............................. 117 Table 5.5 Simulation completion and network latency comparison for a one-hour clock time: MSU HPC Vs. Argonne THETA ............................................................................... 130 Table 6.1 List of the piston skirt measurement locations ........................................................... 143 Table 6.2 Measurement Protocol for the Optical Profilometry .................................................. 144 Table 6.3 Image Analysis Protocol for Omnisurf3D .................................................................. 146 Table 6.4 Comparison of surface properties at side locations for stock (C) and coated (A) piston skirts ........................................................................................................................... 148 ix Table 6.5 Comparing Center Locations for Stock (D) and Coated (B) Piston Skirts ................. 150 Table 6.6 Engine operating conditions ....................................................................................... 153 Table 6.7 Materials of different components .............................................................................. 155 Table 6.8 Mesh details ................................................................................................................ 159 Table 6.9 Performance comparison between worn coated and worn stock conditions .............. 170 Table 7.1 Simulation time and grid details of the cases ............................................................. 179 Table 7.2 Mesh size and no of cells projection for Acadia engine ............................................. 181 x LIST OF FIGURES Figure 1.1 Overview of (a) Power cylinder system (b) 2-D section cut view of Power cylinder system (c) cylinder kit nomenclature. ............................................................................ 1 Figure 1.2 Contribution of cylinder kit on mechanical friction distribution in the power cylinder system............................................................................................................................. 3 Figure 2.1 From Tools to Solution (Flowchart) ............................................................................ 10 Figure 2.2 Schematic View of the Geometry................................................................................ 13 Figure 2.3 Schematic View of the Boundary and Initial Conditions ............................................ 22 Figure 2.4 In-cylinder Pressure Profile of the Engine at 2000 rpm .............................................. 24 Figure 2.5 Prediction of Grid Size Details of Four Simulation Cases: (a) 500 µm (b) 250 µm (c) 125 µm (d) 100 µm ...................................................................................................... 27 Figure 2.6 Averaged Pressure Throughout the Cycle for Different Simulations ......................... 29 Figure 2.7 Circumferential Flow in (a) top land (b)2nd land with top ring end gap view (c) 2nd land with 2nd ring end gap view (d) 3rd land ............................................................... 30 Figure 2.8 Mass Flow Rate at the Outflow Boundary .................................................................. 32 Figure 2.9 Oil Mass Fraction at the Outflow Boundary ............................................................... 32 Figure 2.10 Vector trace at the Inflow Boundary Showing Reverse Blow-by ............................. 34 Figure 2.11 Mass Flow Rate at the Inflow Boundary ................................................................... 34 Figure 2.12 Pressure Across the Piston Length: Comparison Between CASE and Three- Dimensional Simulation at 720 CA ............................................................................. 35 Figure 2.13 Initial Oil Mass Fraction in the Computational Domain at 0 CA (a) liner (b) lands (c) rings (d) grooves .......................................................................................................... 38 Figure 2.14 Comparative Scenario in the Computational Domain at 180 CA (a) lands (b) pressure gradient between third land and oil groove (c) liner (d) rings (e) grooves .... 38 Figure 2.15 Oil Mass Fraction across the Domain at 278 CA (a) lands (b) piston skirt (c) pressure gradient through lands (d) rings (e) grooves (f) liner ................................................... 39 Figure 2.16 Oil Mass Fraction across the Piston and Liner at 360 CA (a) lands (b) rings (c) grooves (d) liner ........................................................................................................... 41 xi Figure 2.17 Oil Mass Fraction across the Piston and Liner at 460 CA (a) lands (b) pressure gradient through lands (c) rings (d) grooves ................................................................ 42 Figure 2.18 Oil Mass Fraction across the Domain at 550 CA (a) lands (b) pressure gradient through lands (c) rings (d) grooves .............................................................................. 43 Figure 2.19 Oil Mass Fraction across the Domain at 720 CA (a) lands (b) rings (c) grooves (d) liner .............................................................................................................................. 45 Figure 2.20 (a) A 2D sketch of the geometry showing the regions of the computational domain (b) oil fractions in different ring regions ...................................................................... 45 Figure 2.21 Oil Evaporation in the Lands throughout the Cycle (a) 180 CA (b) 270 CA (c) 360 CA (d) 460 CA (e) 550 CA (f) 720 CA (g) A 2D sketch of the geometry showing the regions of the computational domain ........................................................................... 46 Figure 2.22 Oil Evaporation in the LR1, LR2 and LR3 Boundaries of the Liner throughout the Cycle (a) 180 CA (b) 270 CA (c) 360 CA (d) 460 CA (e) 550 CA (f) 720 CA (g) A 2D sketch of the geometry showing the regions of the computational domain ................. 47 Figure 2.23 (a) A 2D sketch of the geometry showing the regions of the computational domain (b) oil vapor fractions in different ring regions ............................................................ 47 Figure 3.1 Twist on ring (a) Positive twist (b) Negative twist...................................................... 53 Figure 3.2 From tools to the solution (workflow)......................................................................... 55 Figure 3.3 Inter-ring pressure for second ring input to 3D FEA contact model generated by CASE[29] ..................................................................................................................... 57 Figure 3.4 Second ring motion across the cycle input to 3D FEA contact model (a) ring velocity (b) ring acceleration ..................................................................................................... 58 Figure 3.5 Free shape and deformed ring mesh generated using the ring contact model [39,79] (a) three- dimensional view (b) two-dimensional view ..................................................... 60 Figure 3.6 General trend of twist angle vs. circumferential location: (a) Intake stroke from 0 CA to 180 CA (b) Compression stroke from 181 CA to 360 CA (c) Expansion stroke from 361 CA to 540 CA (d) Exhaust stroke from 541 CA to 720 CA ................................. 61 Figure 3.7 Section cut view of the geometry generated using LINCC [40] ................................. 62 Figure 3.8 Generated mesh of the domain in CONVERGE [46] ................................................. 63 Figure 3.9 Schematic View of the Boundary and Initial Conditions [83] .................................... 64 Figure 3.10 In-cylinder Pressure Profile of the Engine at 2000 rpm [83] .................................... 65 Figure 3.11 Top land Stream traces at a selected crank angle (369 CA) (a) Untwisted xii configuration (b) Twisted configuration ...................................................................... 66 Figure 3.12 second land Stream traces at a selected crank angle (369 CA) (a) Untwisted configuration (b) Twisted configuration. ..................................................................... 67 Figure 3.13 Third land Stream traces at a selected crank angle (369 CA) (a) Untwisted configuration (b) Twisted configuration. ..................................................................... 67 Figure 3.14 Mass flow rate at the outflow boundary for twisted configuration ........................... 69 Figure 3.15 Oil mass fraction at the outflow boundary for twisted configuration ....................... 70 Figure 3.16 Mass flow rate at the inflow boundary for twisted configuration ............................. 71 Figure 3.17 Oil mass fraction at the inflow boundary for twisted configuration ......................... 72 Figure 3.18 Oil consumption vs. brake mean effective pressure at 2000 rpm; comparison with literature ....................................................................................................................... 75 Figure 3.19 A 2D sketch of the geometry showing the regions of the computational domain. ... 76 Figure 3.20 Initial oil mass fraction on the lands at 0 CA for both twisted and untwisted configuration ................................................................................................................ 77 Figure 3.21 Oil mass fraction on the lands at 180 CA (a) untwisted configuration (b) twisted configuration. ............................................................................................................... 78 Figure 3.22 Oil mass fraction on the lands at 278 CA (a) untwisted configuration (b) twisted configuration. ............................................................................................................... 78 Figure 3.23 Pressure on the lands at 278 CA (a) untwisted configuration (b) twisted configuration ................................................................................................................ 79 Figure 3.24 Oil mass fraction on the lands at 550 CA (a) untwisted configuration (b) twisted configuration ................................................................................................................ 79 Figure 3.25 Pressure on the lands at 550 CA (a) untwisted configuration (b) twisted configuration ................................................................................................................ 80 Figure 3.26 (a) Liquid oil mass fractions in different ring regions comparison between twisted and untwisted configuration (b) A 2D sketch of the geometry showing the regions of the computational domain. ........................................................................................... 80 Figure 4.1 Ring pack geometry with (a) no leakage in between second ring and the bottom flank of the groove (b) leakage between second ring and the bottom flank of the groove [29] ...................................................................................................................................... 88 Figure 4.2 Comparison between experimental and CASE predicted results at 1800 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure. ........................................................... 92 xiii Figure 4.3 Comparison between experimental and CASE predicted results at 1800 rpm 51% load (a) 2nd land pressure (b) 3rd land pressure. ................................................................... 93 Figure 4.4 Comparison between experimental and CASE predicted results at 1200 rpm 51% load (a) 2nd land pressure (b) 3rd land pressure. ................................................................... 94 Figure 4.5 Comparison between experimental and CASE predicted results at 1200 rpm 0% load (a) 2nd land pressure (b) 3rd land pressure. ................................................................... 94 Figure 4.6 Comparison between experimental and CASE predicted results at 650 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure. ................................................................... 95 Figure 4.7 Comparison between experimental and CASE predicted results at 650 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure. ................................................................... 95 Figure 4.8 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 1800 rpm 100%) at 1800 rpm 51% load (a) 2nd land pressure (b) 3rd land pressure ........................................................................................................... 96 Figure 4.9 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 1800 rpm 100%) at 1200 rpm 51% load (a) 2nd land pressure (b) 3rd land pressure ........................................................................................................... 97 Figure 4.10 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 1800 rpm 100%) at 1200 rpm 0% load (a) 2nd land pressure (b) 3rd land pressure................................................................................................................. 98 Figure 4.11 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 1800 rpm 100%) at 650 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure ........................................................................................................... 98 Figure 4.12 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 650 rpm 51%) at 650 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure................................................................................................................. 99 Figure 5.1 From tools to the solution (workflow)....................................................................... 103 Figure 5.2 Schematic diagram showing three rings with the corresponding ring face profiles .. 103 Figure 5.3 In-cylinder pressure input to CASE model ............................................................... 105 Figure 5.4 In-cylinder temperature input to CASE model .......................................................... 105 Figure 5.5 Liner temperature ...................................................................................................... 106 Figure 5.6 Cylinder liner deformation ........................................................................................ 107 Figure 5.7. Comparison between experimental and CASE predicted results at 1800 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure .......................................................... 108 xiv Figure 5.8 Temperature boundary conditions of the second ring ............................................... 109 Figure 5.9 Ring cross section temperature distribution .............................................................. 109 Figure 5.10 Second ring motion across the cycle input to 3D FEA contact model (a) ring velocity (b) ring acceleration ................................................................................................... 110 Figure 5.11 Inter-ring pressure for second ring input to 3D FEA contact model generated by CASE ......................................................................................................................... 111 Figure 5.12 Test rig for ring free shape measurement. (a) Initial Test Rig for Ring Tension Measurement (b) Test Rig Modifications for Free Ring Shape Measurement (c) Ring Free Shape Measurement ........................................................................................... 113 Figure 5.13 Second ring free shape profile ................................................................................. 114 Figure 5.14 General trend of twist angle vs. circumferential location........................................ 116 Figure 5.15 Section cut view of the geometry generated using LINCC ..................................... 116 Figure 5.16 Generated mesh of the domain in CONVERGE [46] ............................................. 118 Figure 5.17 Schematic View of the Boundary and Initial Conditions ........................................ 119 Figure 5.18 Section view of the cylinder-kit with initialized oil mass fraction values .............. 120 Figure 5.19 Stream traces at a selected crank angle (418 CA) (a) top land (b) second land (c) third land .................................................................................................................... 121 Figure 5.20 Mass flow rate at the outflow boundary .................................................................. 123 Figure 5.21 Oil mass fraction at the outflow boundary .............................................................. 123 Figure 5.22 Mass flow rate at the inflow boundary .................................................................... 124 Figure 5.23 Comparison between experimental and 3D model predicted results at 1800 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure ................................................ 125 Figure 5.24 Circumferential pressure variation predicted by 3D model (a) 2nd land pressure at 424 CA (b) 3rd land pressure at 493 CA ..................................................................... 126 Figure 5.25 Oil mass fraction on the lands at (a) 0 CA (b) 180 CA (c) 258 CA (d) 361 CA (e) 418 CA (f) 720 CA (g)A 2D sketch of the geometry showing the regions of the computational domain ................................................................................................ 127 Figure 5.26 Oil vapor mass fraction on the lands at (a) 180 CA (b) 258 CA (c) 361 CA (d) 418 CA(e)525 CA (f) 720 CA (g)A 2D sketch of the geometry showing the regions of the computational domain ................................................................................................ 129 xv Figure 5.27 CAs completed vs. no. of cores during a 1-hour clock time at MSU HPC and Argonne THETA ........................................................................................................ 131 Figure 5.28 Average computed network latency vs. no. of cores during a 1-hour clock time at MSU HPC and Argonne THETA .............................................................................. 132 Figure 5.29 CAs completed no. of cores during a 1-hour clock time at Argonne THETA ........ 133 Figure 5.30 Average computed network latency vs. no. of cores during a 1-hour clock time at Argonne THETA ........................................................................................................ 133 Figure 6.1 Selected surfaces for measurement ........................................................................... 142 Figure 6.2 (a) Worn coated piston (b) Worn stock piston .......................................................... 142 Figure 6.3 NPFlex 3D Optical Profiler (Bruker Corporation) [140] .......................................... 144 Figure 6.4 Orientation of the sample for measurement of piston skirt surface........................... 145 Figure 6.5 The surface topology of the worn piston skirt at side location (a) Coated location1 (A1) (b) Stock location 1 (C1) (c) Coated location 4 (A4) (d) Stock location 4 (C4) (e) Worn coated piston, labeled (f) Worn stock piston, labeled. ..................................... 147 Figure 6.6 Comparing the Valleys at Side Locations(a) Valleys at Worn Coated Location 1 (A1) b) Valleys at Worn Stock Location 1 (C1)................................................................. 149 Figure 6.7 Surface topology of the worn piston skirt at center location (a) Coated location 1 (B1) (b) Stock location1 (D1) (c) Coated location 4 (B4) (d) Stock location 4 (D4) (e) Worn coated piston, labeled (f) Worn stock piston, labeled. ............................................... 151 Figure 6.8 Comparing the valleys at center locations (a) Valleys at worn coated location 1 (B1) b) Valleys at worn stock location 1 (D1) ................................................................... 151 Figure 6.9 CASE model setup workflow .................................................................................... 152 Figure 6.10 Combustion gas pressure ......................................................................................... 155 Figure 6.11 Combustion gas temperature ................................................................................... 156 Figure 6.12 Liner temperature .................................................................................................... 156 Figure 6.13 Piston temperature ................................................................................................... 157 Figure 6.14 Piston skirt profile of worn stock and worn coated pistons. (a) vertical trace (b) circumferential trace................................................................................................... 157 Figure 6.15 Piston mesh in CASE .............................................................................................. 159 Figure 6.16 Surfaces showing cylinder-coated piston skirt and cylinder-stock piston skirt ...... 160 xvi Figure 6.17 Piston skirt asperity density (a) worn stock (b) worn coated. ................................. 162 Figure 6.18 Piston skirt asperity radius (a) worn stock (b) worn coated. ................................... 163 Figure 6.19 Piston skirt surface roughness (a) worn stock (b) worn coated. .............................. 163 Figure 6.20 Contact force for coated vs stock piston skirt. ........................................................ 165 Figure 6.21 Hydrodynamic force comparison between the worn coated & worn stock piston skirts. .......................................................................................................................... 165 Figure 6.22 Total force comparison between the worn coated & worn stock piston skirts. ....... 166 Figure 6.23 Minimum skirt oil film thickness comparison between the worn coated & worn stock piston skirts (a) At major thrust side (b) At minor thrust side. .................................. 167 Figure 6.24 Hydrodynamic shear force comparison between the worn coated & worn stock piston skirts. ............................................................................................................... 168 Figure 6.25 Contact friction force comparison between the worn coated & worn stock piston skirts. .......................................................................................................................... 169 Figure 6.26 Total friction force comparison between the worn coated & worn stock piston skirts. .................................................................................................................................... 169 Figure 6.27 Friction power loss comparison between the worn coated & worn stock piston skirt .................................................................................................................................... 171 Figure 6.28 Friction power loss using friction coefficient 0.05 for both stock and coated piston model .......................................................................................................................... 171 Figure 6.29 Friction Power Loss for different projected piston skirt profiles ............................ 172 Figure 6.30 Piston secondary motion comparison between the worn coated & worn stock piston skirts (a)piston tilt (b) piston eccentricity (translation) perpendicular to the wrist pin .................................................................................................................................... 173 Figure 7.1 The three rings oil film thickness prediction by CASE ............................................. 180 Figure 7.2 Difference of the twist angle circumferential variation between the two different ring parts at a selected CA ................................................................................................. 184 Figure 7.3 A flowchart showing cylinder-kit design optimization using 3D multiphysics methodology............................................................................................................... 186 Figure A-1. (a) 3D contour plot at side location 1 (A1) of the coated post run worn piston skirt (b) Zoomed in view of the contour plot (c) Location A1 labeled on the coated post run worn piston skirt......................................................................................................... 189 xvii Figure A-2. (a) 3D contour plot at center location 1 (B1) of the coated post run worn piston skirt (b) Zoomed in view of the contour plot (c) Location B1 labeled on the coated post run worn piston skirt......................................................................................................... 190 Figure A-3. (a) 3D contour plot at side location 1 (C1) of the coated post run worn piston skirt (b) Zoomed in view of the contour plot (c) Location C1 labeled on the coated post run worn piston skirt......................................................................................................... 190 Figure A-4. (a) 3D contour plot at center location 1 (D1) of the coated post run worn piston skirt (b) Zoomed in view of the contour plot (c) Location D1 labeled on the coated post run worn piston skirt......................................................................................................... 191 Figure A-5.GT map of the Cummins R2.8L turbocharged four-cylinder diesel engine............. 192 xviii Chapter 1 INTRODUCTION 1.1 Background and Motivation As increasingly stringent emission standards trigger a switch to a higher number of electrified powertrain architectures in the automotive industry, it is still crucial to continue developing highly efficient, low-emission internal combustion engines to compensate for the limitations of electrical components such as battery energy density or the charging schedule [1]. Cylinder kit tribology has been a major focus in developing engines with reduced friction, lower emission, and higher efficiency. The cylinder kit assembly which is the study subject of our work is a part of the engine power cylinder system as shown in Figure 1.1. Figure 1.1 Overview of (a) Power cylinder system (b) 2-D section cut view of Power cylinder system (c) cylinder kit nomenclature. 1 The internal combustion engine converts thermal energy of the combustible fuel into mechanical energy that moves the piston and eventually the crankshaft. This energy conversion process occurs within the engine power cylinder system. The power cylinder system comprises the following components: piston, piston rings, cylinder liner, wrist pin and connecting rod, see Figure 1.1(a) for three-dimensional view and Figure 1.1(b) for two-dimensional section cut view. In our cylinder kit assembly refers to piston, three piston rings and cylinder liner, as shown in Figure 1.1 (c). The main function of the piston is to transmit the mechanical work done by the ignited gas fuel mixture on the combustion chamber to the crankshaft. The piston presents three grooves on its side to insert the piston rings and the piston outer surfaces between the chamfers are called piston land. The piston skirt is located below the grooves. This surface is in contact with the cylinder liner and its function is to maintain the lateral forces and guide the piston in the cylinder. The three rings of the ring pack are top compression ring, second compression ring and oil ring. The top compression ring is the main component sealing the combustion chamber in order to prevent the high-pressure gas from leaking into the crankcase. The second ring, also called scraper ring, has a hybrid function of scraping the lubrication oil down and sealing the combustion chamber. The oil control ring (OCR) controls the amount of lubrication oil that is supplied to the ring pack. The cylinder of a reciprocating engine is the part through which the piston travels. The cylinder may be sleeved or sleeveless depending on the metal used for the engine block. The sleeves are known as cylinder liners. The nomenclature of the cylinder kit is presented in Figure1.1(c). The cylinder-kit assembly performs two essential functions- sealing and lubrication. It seals the combustion gas chamber from the crankcase by ring motion both in the groove and relative to the piston, and by sliding motion along the cylinder liner. The interaction among liner, piston, and rings during engine operation causes friction and wear between the surfaces. Piston ring friction 2 contributes as much as 25% of total mechanical loss in internal combustion engines. Richardson et al. [2] reported that, in a typical diesel engine about 4 – 15% of the total engine power is wasted as mechanical friction loss. About half of the mechanical friction loss is attributed to the friction in the power cylinder system and the cylinder kit assembly aka the piston and ring pack account for the majority of the friction loss. The detailed breakdown of friction distribution [2] is illustrated in Figure 1.2. Figure 1.2 Contribution of cylinder kit on mechanical friction distribution in the power cylinder system Engine oil consumption is recognized to be a significant source of pollutant emissions in automotive engines. Oil in the exhaust gases contributes directly to particulate emissions. Ingwersen et al. [3] observed that particulate emission from a gasoline-powered engine is increased if oil consumption increases. Trier et al. [4] found that lubricating oil contribution to diesel exhaust emission represented 28% of the total particulates. Sakurai et al. [5] found that at least 95% of the volatile components of nanoparticles and larger particles emitted from a modern, heavy-duty diesel engine under laboratory conditions were unburned lubricating oil. The tradeoff between lowering friction by enabling sufficient oil supply and at the same time minimizing gas leakage and oil consumption calls for a physical understanding of oil transport 3 mechanism in the cylinder-kit assembly. Yet it is difficult and extremely costly to conduct testing and experiments on every series of engines to understand the underlying physics for optimum cylinder kit design. This gives rise to the increasing interest in the development of modeling and simulation tools for better understanding physical processes and providing practical design recommendations for better engine fuel economy, higher efficiency, performance, and endurance. Many researchers have developed numerical models to describe the ring dynamics, piston dynamics, gas flow dynamics and lubrication processes that occur in the assembly. They have analyzed oil transport, oil consumption, blow-by, and underlying mechanisms through numerical models. Ruddy et al. [6] have modeled the hydrodynamic oil film between the ring and liner in an attempt to understand oil consumption. R.J. Gamble et al. [5] developed a model for oil consumption in the ring pack by considering the volume of oil in the piston ring assembly and its residence time in the high-temperature environment. This also considered additional oil transport mechanisms at the ring/cylinder interface and in the piston assembly, as well as the worn condition of the components and the shape of the oil film accumulating on the cylinder wall. Tian et al. [7– 9] developed ring pack models that predict ring dynamics and their effect on blow-by and oil transport along the piston and liner. Thirouard et al. [13–15] presented a global oil transport scheme where he investigated different oil transport mechanisms across the piston. Yilmaz [13] pointed out important oil consumption characteristics and sources in SI engines for different engine operating conditions via a multi-species liner evaporation model. Baelden et al. [14] used a multi- scale curved beam finite element model to explore oil transport at the ring-liner and ring-groove interfaces. This model was extended by Liu [15] to find ring structural response to ring-liner as well as ring-groove interaction. Later, Bhouri [16] adapted the model to study the impact of thermal distortion and expansion on ring-bore conformability and twist calculation. Wang [17] 4 studied three different flow regimes of blow-by and their effect on oil transport. Fang [18] extended the quantitative models of oil transport mechanism with an additional multi-phase CFD (Computational Fluid Dynamics) model to investigate oil transport from piston land to liner. Such transport was defined as “bridging” in [19]. There have only been a few studies about multi-dimensional, multi-phase CFD simulation of the cylinder kit in the past. Felter [20] developed a compressible Navier-Stokes equation-based 2D (two- dimensional) free surface model which focuses on calculating oil film outside the piston ring. The domain of the model was only around the ring. Shahmohamadi et al. [21,22] introduced an isothermal, two-phase, two-dimensional CFD model where the computational domain is limited to the compression ring running face, not the entire ring pack. The key finding of the study was the impact of ring running face on cavitation and other tribological conditions. Hronza and Bell [23] introduced a 2D, multi-phase CFD model using Navier–Stokes equations instead of the Reynolds equation to model the oil phase. The model considered hydrodynamic lubrication, oil transport, and radial ring dynamics, including the effects of surface tension and wetting angle, to calculate oil distribution in leading and trailing areas surrounding the piston ring. Veettil and Shi [24] developed a two-phase, three-dimensional CFD model where the computational domain is the piston from the top land to the piston skirt and the liner. They studied the effect of the position of the rings in the grooves and the position of the end gaps on blow-by and oil consumption. Due to the complexity of ring dynamics, this was not considered in the model. Olivia et. al. [25] presented a two-dimensional CFD model to calculate the gas flow through the piston ring pack. They used a 1D (one-dimensional) tool to adjust their CFD method by adjusting the leakage of each ring with respect to its groove sides. The influence of the ring end gaps was also compensated with adjusted axial ring motion. They further extended this effort by adding oil to explore oil transport 5 mechanisms throughout the cycle. The effect of different oil filling ratios in the oil ring reservoir was also studied [26]. Bartel et al. [27] developed a multi-phase, 2D CFD simulation model considering piston ring dynamics, solid contact, mixed friction, and transient boundary conditions for combustion chamber pressure and temperature and thermal distortion on piston and liner. Urzua et al. [28] conducted a 3D CFD simulation of the area between the top of combustion chamber and the first ring groove during compression stroke. The focus of the study was the fuel/oil entry in the combustion chamber to explore the effect of top ring up scraping on low speed pre-ignition (LSPI). The Energy & Automotive Research Laboratory at Michigan State University has been contributing to cylinder kit modeling in collaboration with Mid – Michigan Research. These groups developed a quasi-one-dimensional system of programs, CASE (Cylinder-kit Analysis System for Engines) [29]. Brombolich [30–32] developed a ring pack analysis methodology based on a concept of a quasi-steady orifice-volume gas flow model. Ejakov et.al. [33] developed a three- dimensional beam model of a piston ring. This model established an interface between the in- cylinder combustion model and the piston ring. They simulated the ring pack kinematics and gas dynamics for a deactivated cylinder [34]. They also modeled the three-dimensional ring twist in the piston grooves using space beam elements [35]. Chui et al. [36] proposed a numerical methodology to develop a three-dimensional model of the top compression ring lubrication that accounts for the partially flooded condition, coupled with the oil evaporative effect. They also developed a three-dimensional piston elastohydrodynamic (EHL)lubrication model, using the finite element approach to simulate the dynamics of a piston that allows its structure to deform elastically within the cylinder bore over an engine cycle [37]. Panayi et al. [38] developed a novel, three-dimensional model, considering the secondary motion in the wrist-pin plane to predict piston 6 dynamics in addition to the axial and thrust plane motions of the piston. Cheng et. al. [39] developed a three-dimensional piston ring model, using a finite element method that predicts the piston ring conformability with the cylinder wall separation gap between the interfaces and the interaction between the ring and piston groove sides. Kharazmi [40] developed an automated cylinder-kit geometry generating program known as LINCC [Link CASE to Computational Fluid Dynamics (CFD)]. This considered the complicated geometrical details of the ring pack such as thermal distortion of piston and liner, ring twist and ring/groove conformability. He also developed a three-dimensional CFD model to analyze the gas flow between the cylinder liner and the piston. The challenges of developing a multi-phase three-dimensional model are manifold. Multiple physical processes of gas and oil transport take place simultaneously or concurrently requiring very fine meshing, small timestep, and tight convergence tolerance to resolve these phenomena; hence making the problem complex, time-consuming and resource constrained. With the enormous progress in the field of high-performance computing in the last few decades, there has been a huge improvement in processing time and capability, so that models are more effectively applied to today’s decision-making and design. This cuts down the lead time and costs for the assessment of cylinder kit performance and optimization of designs more than ever. The current work is aimed at establishing a multi-phase, three- dimensional, CFD analysis methodology to study the blow-by and oil consumption characteristics of the cylinder kit throughout the cycle using computational tools and high-performance computing. This is the very first attempt of demonstrating gas-oil transport in the cylinder kit using a three-dimensional model. This model proposes the oil transport and oil evaporation mechanisms inside the cylinder kit assembly. The oil supply mechanism to the cylinder kit assembly is also investigated. Later, the model is used to explore the effect of piston ring twist and leakage area on oil-gas transport. 7 Coating materials have been gaining interest among the piston manufacturers to minimize wear and friction in the engine. In the later part of the present work, a novel abradable powder coating surface is characterized and compared with stock coatings. Finally, the performance of piston skirts with the two coatings are assessed using CASE. 1.2 Structure of Dissertation The dissertation is organized as follows- In chapter 2, the development methodology of the three-dimensional multiphase model has been demonstrated. A grid sensitivity test is carried out. Blowby, reverse blowby are estimated and compared with CASE results. Oil and gas transport are discussed. Chapter 3 presents the development work regarding the dynamically twisted second ring analysis. A three-dimensional model is developed for the twisted geometries across the cycle.The twisted and untwisted configurations are compared. In chapter 4, an optimization study is presented using CASE and HEEDS on a Cummins Acadia engine. Chapter 5 illustrates the experimental validation of the three-dimensional methodology on a Cummins Acadia engine. In Chapter 6, the tribological studies on a novel abradable powder coating is presented. Chapter 8 provides the concluding remarks regarding the modeling effort and recommends steps for future work. 8 Chapter 2 THREE-DIMENSIONAL MULTI-PHASE PHYSICS-BASED MODELING METHODOLOGY TO STUDY ENGINE CYLINDER-KIT ASSEMBLY TRIBOLOGY AND DESIGN CONSIDERATIONS- PART I 2.1 Introduction Understanding cylinder-kit tribology is pivotal to durability, emission management, reduced oil consumption, and efficiency of the internal combustion engine. This work addresses the understanding of the fundamental aspects of oil transport and combustion gas flow in the cylinder kit, using simulation tools and high-performance computing. A dynamic three-dimensional multi- phase, multi-component modeling methodology is demonstrated to study cylinder-kit assembly tribology during the four-stroke cycle of a piston engine. The percentage of oil and gas transported through different regions of the piston ring pack is predicted, and the mechanisms behind this transport are analyzed. The velocity field shows substantial circumferential flow in the piston ring pack, leading to blowback into the combustion chamber during the expansion stroke. Oil initialization and management of a continuous supply of oil throughout the cycle are observed to govern how much oil would be lost to the crankcase and combustion chamber. The calculated blow-by results agree with the results of a quasi-one-dimensional cylinder-kit analysis system of programs known as CASE (Cylinder-kit Analysis System for Engines). Implementing this three- dimensional methodology leads to a better understanding of cylinder-kit fluid flow physics. The findings presented in this work pave the way to further the ongoing development effort of optimum cylinder kit designs with controlled gas leakage, low oil consumption, and low cylinder kit friction. 9 2.2 Workflow As a first step, a quasi-one-dimensional model is developed in CASE. The results are then used to generate the geometry using LINCC. Finally, the results from CASE and the geometry from LINCC are used to develop a multi-phase three-dimensional model in CONVERGE. The detailed workflow is shown in Figure 2.1. A grid sensitivity test is performed and from the comparison of the grid resolution, the size of one grid is selected for further analysis. Blow-by and pressure distribution along the piston length at the end of the cycle found from the simulation are compared to the CASE results. The physical phenomena across the cycle, namely blow-by, reverse blow-by, circumferential flow across the lands, liquid oil, and oil vapor distribution across the regions throughout the cycle are finally discussed. Figure 2.1 From Tools to Solution (Flowchart) 10 2.2.1 One-Dimensional Simulation & Geometry Generation The first step in the workflow is developing a one-dimensional simulation model. Therefore, the Cylinder-kit Assembly System for Engines (CASE) is used to model the cylinder kit components [29]. CASE includes the lubrication model using Reynolds’ equation [41], hydrodynamic contact model using Patir-Cheng’s average flow model [42], asperity contact model using Greenwood & Tripp equation [43], Wear model using Archard’s Model [44], and the oil consumption model using Antoine’s equation [45]. The competing forces acting on the ring axially and radially such as gas pressure force, inertia force, friction force, lubricant squeeze force and lubricant adhesive force, ring tension and the competing forces acting on the piston such as inertia, hydrodynamics, contact, shear, frictional forces are considered in CASE. Note that, for lubrication CASE assumes a fully flooded condition for all three rings. The one-dimensional analysis provides the necessary boundary conditions such as crank-angle- resolved motion data for piston, ring, ring-bore conformability, pressure, etc. for the subsequent three-dimensional analysis. The cylinder-kit assembly in this study is comprised of two compression rings: a keystone top ring, a tapered front face second ring, an oil control ring and a simplified cylindrical piston skirt in a small-bore experimental engine, as shown in the two- dimensional sketch in Figure 2.2. Some of the key input parameters for the CASE model include geometrical dimensions such as detailed measurements of the piston, rings, and liner, including end clearance of piston ring pack, location of ring end gaps; material characterization data such as material properties, surface characteristics, wear & friction coefficients; thermodynamic attributes, for example, combustion chamber pressure curve and temperature of the piston, liner, and intake combustion gas temperature; flow characteristics, for example, fluid properties of the engine oil, flow coefficients, 11 etc. Some of the engine operating conditions are presented in Table 2.1, as follows. Table 2.1 Engine operating condition Parameter Value Bore 50.6mm Stroke 50mm Connecting rod length 130 mm Engine speed 2000 rpm Compression ratio 14 Engine oil SAE 10W-30 Based on these inputs, a simulation of the cylinder kit is performed, and the results are used as boundary conditions for the next steps, i.e., the geometry generation and three-dimensional model development. The next step toward three-dimensional simulation is geometry generation. A custom program, LINCC (Link CFD CASE) developed in [40], is used to establish a bridge between the CASE and the 3D CFD model. LINCC creates the ring pack geometry in STL (Stereolithographic) format progressively. At first, it creates the rings, the piston, and the cylinder liner. The construction of the separate geometries is based on CASE input parameters and crank-angle-resolved ring location, considering ring-bore conformability as well as the separation gap between the interfaces as described in [39]. Next, they are assembled to form the geometry of the current study. Since this is the first attempt to develop a multi-phase, three-dimensional model, the deformation of the piston or ring twist is not considered in this study to avoid its additional complexity in physics. The results from the CASE simulation and the geometry from LINCC are then used for three- 12 dimensional cylinder kit multi-phase model development, using CONVERGE [46]. 2.3 Three-Dimensional Multi-phase Model Development Figure 2.2 Schematic View of the Geometry The computational domain starts at the top land of the top compression ring and ends at the piston skirt. The combustion chamber and the crankcase are not included in the domain. A two- dimensional sketch and a three-dimensional half-cut view of the geometry are illustrated in Figure 2.2. In this section, the sub-models and governing equations used to develop the three- dimensional multi-phase model are explained. The boundary and initial conditions are also described. 2.3.1 Volume of Fluid (VOF) Model The mixture model [47] with the volume of the fluid (VOF) modeling approach [48] is used for this analysis. As the name implies, the VOF method tracks the volume of fluid within each cell 13 [49]. The volume of fluid is represented by the void fraction 𝛼, which is the fraction of the cell's volume that does not contain fluid: 𝛼=0 (The cell contains only liquid), 0 < 𝛼 < 1 (The cell contains both liquid and gas) , (2.1) 𝛼=1 (The cell contains only gas). This value is tracked throughout the domain. In the present study, the gas is air and the liquid is oil. SAE 10-W30 is chosen as oil which is assumed to be incompressible, and the compressible gas is assumed to be a combination of only 𝑁2 and 𝑂2. It is assumed that 𝐶7 𝐻16 is the vapor constituent of liquid oil, which will be discussed later in this section. The local value of the void fraction does not contain any information about the shape or location of any interface within the cell. In the current study, the High-Resolution Interface Capturing (HRIC) scheme [50] is used to construct an interpolated curved interface between air and oil. This scheme blends up-winding and down-winding trends to balance stability and accuracy. It maintains a sharper interface than a VOF simulation without an interface capturing scheme. The blending factor is determined by numerical parameters. However, it is computationally expensive. The VOF approach follows the Individual Species Solution Method [49] for compressible air and incompressible engine oil. In this method, CONVERGE transports species, momentum, and energy, and calculates the void fraction based on the species mass fraction. 14 The void fraction is solved with the following conservation equation: 𝜕𝛼 𝜕(𝛼𝑢𝑖 ) + =0 (2.2) 𝜕𝑡 𝜕𝑥𝑖 Global density (𝜌) is computed as 𝜌 = 𝛼𝜌𝑔 + (1 − 𝛼)𝜌𝑙 (2.3) where 𝜌𝑔 and 𝜌𝑙 stand for gas density and liquid density in the cell, respectively. The void fraction is not transported directly, as in Equation (2.2). At first, the species is solved using the following species transport equation: 𝜕𝜌𝑚 𝜕𝜌𝑚 𝑢𝑖 𝜕 𝜕𝑌𝑚 + = (𝜌𝐷 ) , 𝑚 = 1, … 𝑛 (2.4) 𝜕𝑡 𝜕𝑥𝑖 𝜕𝑥𝑖 𝜕𝑥𝑖 where 𝜌𝑚 is the density of species m, 𝑌𝑚 is the mass fraction of species m, 𝐷 is the mass diffusion coefficient, 𝑢𝑖 is the velocity component in the ith direction, and 𝑛 is the total number of species. In the current study, there are four species: one liquid species (oil), and three gaseous species (𝑁2 ,𝑂2 and C7 H16 ). Then the void fraction is computed, using the species mass fractions in the cell. The solver sums the mass of gas species in the cell to calculate the total gas mass in a cell. The remainder of the mass in a cell is considered the liquid mass. In other words, 15 𝑛𝑔 𝑚𝑔 = ∑ 𝑌𝑚 𝑚=1 (2.5) 𝑚𝑙 = 1 − 𝑚𝑔 where 𝑚𝑔 is the total gas mass fraction, 𝑚𝑙 is the total liquid mass fraction, and 𝑛𝑔 is the total number of gas species (three in the present work). Finally, the void fraction is expressed as 𝑚𝑔 𝜌𝑔 𝛼= 𝑚𝑔 𝑚𝑙 (2.6) ( + ) 𝜌𝑔 𝜌𝑙 2.3.2 Governing Equations & Solver Setup The governing equations are as follows. The mass equation is 𝜕𝜌 𝜕𝜌𝑢𝑖 + = 𝑆𝜌 (2.7) 𝜕𝑡 𝜕𝑥𝑖 The momentum equation is 𝜕𝜌𝑢𝑖 𝜕𝜌𝑢𝑖 𝑢𝑗 𝜕𝑃 𝜕𝜎𝑖𝑗 + = − 𝜕𝑥 + + 𝑆𝑚𝑜𝑚 (2.8) 𝜕𝑡 𝜕𝑥𝑗 𝑖 𝜕𝑥𝑗 16 The energy equation is 𝜕𝜌𝑒 𝜕𝜌𝑢𝑖 𝑒 𝜕𝑢 𝜕𝑢 𝜕 𝜕𝑇 + = −𝑃 𝜕𝑥𝑖 + 𝜎𝑖𝑗 𝜕𝑥 𝑖 + 𝜕𝑥 (𝐾 𝜕𝑥 ) 𝜕𝑡 𝜕𝑥𝑖 𝑖 𝑗 𝑖 𝑖 𝜕 𝜕𝑌𝑚 + 𝜕𝑥 (𝜌 ∑𝑚 𝐷ℎ𝑚 ) + 𝑆𝑒 (2.9) 𝑖 𝜕𝑥𝑖 The viscous stress tensor (𝜎𝑖𝑗 ) in the equations above is given by 𝜕𝑢 𝜕𝑢 2 𝜕𝑢 𝜎𝑖𝑗 = 𝜇 (𝜕𝑥 𝑖 + 𝜕𝑥𝑗) + (𝜇 ′ − 3 𝜇) (𝜕𝑥𝑘 𝛿𝑖𝑗 ) (2.10) 𝑗 𝑖 𝑘 where 𝜌𝑚 is the density of species m, 𝜎𝑖𝑗 is the stress tensor, 𝑌𝑚 is the mass fraction of species m, 𝐷 is the mass diffusion coefficient, ℎ𝑚 is the enthalpy of species m, 𝑒 is specific internal energy, 𝑆 is the source, 𝑢𝑖 is the velocity component in the ith direction, μ' is the dilatational viscosity (set to zero), and 𝛿𝑖𝑗 is the Kronecker delta. The Navier-Stokes solver [51][49] for the volume of the fluid model is selected. A Reynolds- averaged Navier-Stokes (RANS) turbulence model of standard k-ε [52] is utilized to model the turbulence of the flow. The dissipation that is introduced in the flow equations by the current turbulence model helps the convergence of the results by removing the low-frequency fluctuations of the simulation runs. Since the RANS-standard k-ε turbulence model is activated in the study, the molecular viscosity is augmented by a turbulent viscosity term. The RANS model represents the turbulent viscosity (𝜇𝑡 ) as a function of k and ε 𝑘2 𝜇𝑡 = 𝑐𝜇 𝜌 (2.11) 𝜀 17 where 𝑐𝜇 is a turbulence model constant, 𝑘 is the turbulent kinetic energy, and 𝜀 is the turbulent dissipation. Again, due to the turbulence model, the conductivity (K) in equation (2.9) is replaced by turbulent conductivity given by: 𝜇𝑡 𝐾𝑡 = 𝐾 + 𝑐𝑃 (2.12) 𝑃𝑟𝑡 Here, 𝑃𝑟𝑡 is turbulent Prandtl number and 𝜇𝑡 is the turbulent viscosity. The Redlich-Kwong cubic equation of state [53] for the species is employed to couple the density, pressure and temperature equations. The general equation for multi-species simulation can be expressed as 𝑅𝑇 𝑎𝑚 𝑃= − 2 2) (2.13) 𝑣 − 𝑏 (𝑣 + 𝑢𝑏𝑚 𝑣 + 𝑤𝑏𝑚 where P is pressure, T is temperature, R is constant and the rest of the terms are coefficient, details of which can be found in [49]. The pressure is not calculated directly from the equation of state but instead is used indirectly in the PISO (Pressure Implicit with Splitting of Operator) algorithm [54] to ensure that the equation of state is satisfied. This algorithm uses the elliptic nature of the pressure equation to transmit information more quickly through the domain. This is a predictor-corrector method of solving pressure-velocity coupling. PISO significantly accelerates the convergence of the compressible solver. The PISO algorithm as implemented in CONVERGE starts with a predictor step where the momentum equation is solved. After the predictor, a pressure equation is derived and solved, which 18 leads to a correction, which is applied to the momentum equation. This process of correcting the momentum equation and re-solving can be repeated as many times as necessary to achieve the desired accuracy. After the momentum predictor and first corrector step have been completed, the other transport equations are solved in series. It is necessary to solve the mass and momentum equations together for the proper calculation of the pressure gradient in the momentum equation [49]. To balance simulation speed, accuracy, and stability, convergence tolerance values presented in Table 2.2 are used while solving the equations. Table 2.2 Convergence tolerance values Parameter Convergence Tolerance Momentum 1.0e-6 Pressure 1.0e-8 Density 1.0e-5 Energy 1.0e-5 PISO 1.0e-02 PISO convergence criterion multiplier 100 The transient solver with full hydrodynamic simulation mode is chosen. The simulation timestep is set to vary between 1e-08 seconds and 0.0001 seconds, which is adjusted using a variable time step algorithm based on CFL (Courant–Friedrichs–Lewy) numbers [55]. The CFL numbers estimate the number of cells through which the related quantity will move in a single time-step. A higher number generally yields a less computationally expensive simulation. The convective CFL number (𝑐𝑓𝑙𝑢 ), the speed of sound CFL number (𝑐𝑓𝑙𝑚𝑎𝑐ℎ ), and the diffusive CFL number (𝑐𝑓𝑙𝜈 ) are given respectively as: ∆𝑡 𝑐𝑓𝑙𝑢 = 𝑢 (2.14) ∆𝑥 19 ∆𝑡 𝑐𝑓𝑙𝑚𝑎𝑐ℎ = 𝑐 (2.15) ∆𝑥 𝛥𝑡 𝑐𝑓𝑙𝜈 = 𝜈 (2.16) 𝛥𝑥 where ∆t is the time-step, ∆x is the grid spacing, u is the cell velocity, c is the speed of sound, and 𝜈 is the viscosity. The equations are then solved by marching in time and the results are calculated for each timestep. This allows resolving the oil-air distribution throughout the domain in 3D at every crank angle. However, it is very time-consuming to resolve the whole domain with a transient solver. To speed up the solving process a little bit, the CFL numbers are set to the upper limit for the volume of fluid approach, namely, maximum convection CFL limit 1, maximum diffusion CFL limit 2, maximum Mach CFL limit 10. 2.3.3 Homogeneous Relaxation Model A homogeneous relaxation model is used to analyze the evaporation of liquid oil. The multi- species nature of the oils increases the complexity of the evaporation process of an oil film. Because hydrocarbon species have different volatilities, their rate of evaporation from the film surface differs considerably [56]. The transport of each species through the film to the surface is governed by transient diffusion. Depending on the application, this transport could result in a non- uniform distribution of the hydrocarbon species within the liquid film. Moreover, the composition of the film is unsteady and changes with the liner location. This characteristic of the oil film directly influences the oil evaporation since the vapor mixture in the gas side depends on the liquid oil composition at the interface. Since the exact composition of oil is unavailable, it is assumed 20 that the most volatile constituent of the oil is n-heptane (C7 H16 ). This is considered to be the liquid- vapor (gas species) of liquid oil SAE 10W-30 in the evaporation model known as the homogeneous relaxation model [57][58]. The model describes the rate at which the instantaneous mass approaches its equilibrium value. The governing equations are as follows. 𝑑𝑥 𝑥̅ − 𝑥 = (2.17) 𝑑𝑡 𝜃 𝜃𝐸 = 𝜃0 𝛼 −0.54 𝜑 −1.76 (2.18) 𝑃𝑠𝑎𝑡 − 𝑃 𝜑= (2.19) 𝑃𝑐 − 𝑃𝑠𝑎𝑡 Here, 𝑥 is species mass fraction, 𝑥̅ is equilibrium mass fraction, 𝜃 is time-scale over which x relaxes to 𝑥̅ , 𝜃𝐸 is time-scale for evaporation, 𝜃0 is time scale coefficient, 𝛼 is the void fraction, 𝜑 is non-dimensional pressure ratio, 𝑃𝑠𝑎𝑡 is saturation pressure, 𝑃𝑐 is critical pressure, and 𝑃 is pressure. 21 2.3.4 Boundary & Initial Conditions Figure 2.3 Schematic View of the Boundary and Initial Conditions In CONVERGE each boundary, regardless of type, is assigned a boundary condition for the partial differential conservation equations as either Dirichlet (specified value) or Neumann (specified first derivative value) [49]. Mathematically, these boundary conditions are given as: 22 𝜙 = 𝑓(𝐷𝑖𝑟𝑖𝑐ℎ𝑙𝑒𝑡) and (2.20) 𝜕𝜙 = 𝑓(𝑁𝑒𝑢𝑚𝑎𝑛𝑛) 𝜕𝑥 where 𝜙 is a solved quantity (e.g., pressure, energy, velocity, or species) and 𝑓 is the specified value or specified derivative on the boundary. A schematic view of the boundary conditions used in this problem is illustrated in Figure 2.3. For inflow and outflow boundary conditions, pressure and temperature from CASE are applied at the inlet and outlet of the cylinder-kit assembly as shown in Figure 2.3. The in-cylinder pressure profile throughout the cycle calculated using CASE is presented in Figure 4. The Dirichlet pressure and temperature boundary conditions applied here are total pressure (𝑃𝑡𝑜𝑡𝑎𝑙 ) and total temperatures (𝑇𝑡𝑜𝑡𝑎𝑙 ) respectively. The corresponding static pressure is calculated from the following relation: −1 𝛾 − 1 𝑢𝑖2 𝑃𝑠𝑡𝑎𝑡𝑖𝑐 = 𝑃𝑡𝑜𝑡𝑎𝑙 (1.0 + ) (2.21) 2 𝛾𝑅𝑇 And, the corresponding static temperature is determined from: −1 𝛾 − 1 𝑢𝑖2 𝑇𝑠𝑡𝑎𝑡𝑖𝑐 = 𝑇𝑡𝑜𝑡𝑎𝑙 (1.0 + ) (2.22) 2 𝛾𝑅𝑇 where 𝑃𝑠𝑡𝑎𝑡𝑖𝑐, 𝑇𝑠𝑡𝑎𝑡𝑖𝑐 , 𝛾 and R are static pressure, static temperature, the ratio of specific heat, and ideal gas constant, respectively. 23 Figure 2.4 In-cylinder Pressure Profile of the Engine at 2000 rpm Turbulent kinetic energy (TKE) and turbulent dissipation (eps) are set as default Dirichlet boundary conditions for the inflow and outflow boundaries. The fixed surface movement and stationary wall motion boundary conditions are applied at cylinder liner, piston, and ring surfaces where velocity is assumed to be zero, hence no-slip wall boundary condition. Since the k-ε turbulence model is used, the law of the wall boundary condition (Launder-Spalding model) is used [59]. Since the near-wall resolution is not adequate in this study, the temperature law of the wall [52] boundary condition is used. For all boundaries, the temperature boundary condition at different boundaries found from CASE is presented in Table 3. The O'Rourke and Amsden model [60] is used for wall heat transfer. 24 Due to the time and resource constraints associated with carrying out a fine mesh (micron scale) three-dimensional simulation with a moving boundary, the piston-ring movement is not modeled in CONVERGE. However, all the boundary conditions mentioned above are found from CASE simulation which considered this. Table 2.3 Temperature boundary condition Boundary Temperature (K) Inflow gas, 𝑇𝑔𝑎𝑠𝑖𝑛𝑡 330 Lands- Top; Second; Oil 863; 623; 471 Grooves - Top; Second; Oil 793; 491.6; 453 Rings-Top; Second; Oil 828; 557.3; 462 Piston Skirt 403 Liner 403 Outflow,𝑇𝑠𝑢𝑚𝑝 300 Due to calculation time and resource restrictions, it is not viable to calculate multiple cycles of the engine. It is assumed that the results would repeat periodically every cycle. The oil is initialized at the start of the simulation. To initialize oil and gas, the computational domain is divided into four regions, namely: top ring region, second ring region, oil ring region, piston skirt region as shown in Figure 2.3. The boundaries associated with each region are listed in Table 2.4. Table 2.4 Boundaries assigned to the regions Regions Boundaries Top ring region -Top land -Top ring groove -Top ring -LR1(part of the liner adjacent to the top ring) Second ring region -Second land -Second ring groove -Second ring -LR2 (part of the liner adjacent to the second ring) Oil ring region -Third land -Oil ring groove -Oil ring -LR3 (part of the liner adjacent to oil ring) Piston skirt region -Piston skirt -Lskirt (part of the liner adjacent to piston skirt) 25 Since the main oil source is the crankcase, the oil is initialized gradually, increasing from the piston skirt. The surface roughness of piston skirt, liner and all three rings are 10 µm, 2 µm and 0.4 µm, respectively. This is not considered in the assumption, rather a gradually increasing oil mass fraction from skirt region initialization approach is taken. It is assumed that there is 5%,10%, and 15% and 100% of oil mass uniformly initialized in the top ring region, second ring region, oil ring region, and piston skirt region, respectively. Variation of initial oil mass fraction in these regions has not been studied in the present work. The remaining fluid domain is initialized with gas (N2 and O2 ). However, in practice, the product of combustion primarily contains H2 O, CO2 andN2 . The regions, associated boundaries, and the corresponding initial conditions are depicted in the two-dimensional schematic view of the geometry in Figure 2.3. 2.4 Results and Discussion 2.4.1 Grid Sensitivity Test Based upon the numerical setup described in the previous section, a mesh grid sensitivity analysis is carried out to evaluate the convergence of numerical solutions. The performance of the CFD models is determined by the time required by the simulations and the accuracy of the results. These are heavily affected by the choice of the mesh size. Thus, four simulations are carried out. They are based on identical boundary conditions, initial conditions, and properties, but set up using four different mesh sizes. All the simulations in the current study are performed in the High-Performance Computing Center (HPCC) of Michigan State University. Jobs are submitted to the system, and computational resources are provided, based on availability. HPCC is comprised of different clusters. Each cluster runs the Community Enterprise Operating System (CentOS) 7, each with a different combination 26 of memory and interconnects. Each cluster uses the SLURM resource manager for job management. Grid size details, with total simulation time, for the four cases are specified in Table 2.5. In the table, dx, dy, and dz indicate cell dimensions in the x, y, and z directions, respectively. Figure 2.5 Prediction of Grid Size Details of Four Simulation Cases: (a) 500 µm (b) 250 µm (c) 125 µm (d) 100 µm Table 2.5 Simulation time and grid details of the cases Case no. dx(µm) dy(µm) dz(µm) Total no. of cells Total CPU hours 1 500 500 500 25,238 750 2 250 250 250 171,435 74,800 3 125 125 125 1,266,257 443,000 4 100 100 100 2,403,538 522,000 The grids used to generate solutions must be sufficiently fine to resolve the oil-air interface. Moreover, the grids need to be small enough to capture the tiny ring-groove clearances on the micron scale. Since the interface geometry is unknown, all four cases have uniformly distributed 27 grid sizes for each cell. To compare the four cases, the grids are inspected at the same location and crank angle degree. Since it is hard to illustrate and visualize the difference in grid sizes in a three- dimensional picture, 2-D images of the oil ring groove along the z-x co-ordinate system at 720 CA is presented here. Figure 2.5(a), Figure 2.5(b), Figure 2.5(c), and Figure 2.5(d) refer to the oil mass fraction predicted by grid sizes 500 µm, 250 µm, 125 µm, and 100 µm, respectively. It is observed that the coarsest mesh is unable to capture the oil mass fraction, as shown in Figure 2.5(a). It predicts the oil mass fraction to be 0.05. As the mesh is refined more and the number of cells increases, the oil mass fraction is more precisely captured. The 250 µm mesh captures a slight gradation in oil mass fraction on the range of 0.15 to 0.45; refer to Figure 2.5(b). The 125 µm mesh predicts the oil mass fraction to be 0.85, but it is unable to capture any gradation. Among the four grid sizes, the 100 µm demonstrates the best resolution, as can be seen in Figure 2.5(d). The 100 µm captures the oil mass fraction to be 0.85 at the left corner as in the 125 µm case, but in addition to that, it captures the gradation of the oil mass fraction. The averaged pressure along the entire ring pack throughout the cycle for different grid sizes is also plotted in Figure 2.6. As shown in Figure 2.6, it is observed that the refinement of the grid 2 size shows considerable changes (10 𝑀𝑃𝑎 𝑜𝑟 29.0075 𝑝𝑠𝑖) in pressure. Both Figure 2.5 and Figure 2.6 imply that a grid-independent solution demands a very refined grid, which would require more computational resources and time to run. Due to time and resource constraints, the 100 µm case is considered for further investigation. 28 Figure 2.6 Averaged Pressure Throughout the Cycle for Different Simulations 2.4.2 Circumferential Flow The present work resolves the 3D effect of oil-gas flow inside the cylinder kit for the first time. The circumferential flow in the domain of the problem is investigated in different regions of the piston ring pack by computing the velocity field across the piston ring. The impact of the circumferential flow among the rings, grooves, and lands can greatly impact blow-by oil transport, which will be discussed later. To show the general circumferential flow behavior, the circumferential flow on the three lands of the studied cylinder-kit geometry at a selected CA (Crank Angle) 388 is depicted in Figure 2.7. It is observed that the flow velocity across the lands is influenced by the ring end gap. In the top land, flow is mostly axial, and it only becomes partially circumferential closer to the ring end gap which is shown in Figure 2.7(a). In the top land, small vortices ranging from 10 to 60 m/s are formed. Figure 2.7(b) and Figure 2.7(c) show the flow in the second land with a view of the top 29 ring end gap and second ring end gap, respectively. The flow in the second land is greatly affected by the respective location of the first and second ring end gaps. Because of the small clearance between the surface of the second ring and its groove, small vortices ranging from 10-120 m/s form on the 2nd land. As shown in Figure 2.7(d), the same conditions are true for the 3rd land. The vorticity of the flow exiting the second ring end gap into the third land is much larger, but only one vortex forming on each side of the 2nd ring end gap is observed. The above results suggest that the geometry of the piston ring plays an important role in determining the circumferential flow [40]. Figure 2.7 Circumferential Flow in (a) top land (b)2nd land with top ring end gap view (c) 2nd land with 2nd ring end gap view (d) 3rd land 30 2.4.3 Estimating Blow-by Blow-by is the gas flow from the combustion chamber to the crankcase. This flow can occur through leakage passages such as ring end gaps; circumferentially through the lands; and through the gap between the liner and ring front or backside of the ring when the ring loses its stability [24,61–64]. Blow-by is highly critical to engine cylinder kit design. Calculation of blow-by requires measurement of the amount of the gas that passes from the combustion chamber through the ring pack into the crankcase during a complete cycle. Blow-by is calculated by summing the mass flow rate that goes through the outlet boundary over the whole cycle. Figure 2.8 shows the mass flow rate at the outlet boundary throughout the cycle found from the simulation. In this study, the engine RPM is set to 2000 and thus it takes 0.06 seconds to complete a cycle. 1 Therefore, 1 CA corresponds to 12 millisecond. The discharged mass in kg is found by integrating the area under the mass flow rate versus time graph. The discharged mass is calculated to be 0.0028 kg. This cumulative mass is the sum of the mass of oil and mass of gas. Observing the oil mass fraction at the outflow boundary in Figure 2.9 the averaged oil mass fraction is found to be 0.99. Considering the gas mass fraction to be 0.01 at the outlet boundary results in discharged gas per cycle to be 4.79 ∗ 10−5 kg/s. The final step to calculate blow-by is to convert the mass flow rate from kg/s to liter/min. It is kg assumed that air density is 1.1766 m3 , considering air pressure and temperature to be 101325 Pa and 300K, respectively. Finally, the blow-by is calculated to be 2.477 liters/min. The one- dimensional simulation in CASE estimates average blow-by through oil ring to be 2.81 liters/min. Assuming blow-by through the oil ring to be equal to blow-by through the outflow boundary, it can be said that the three-dimensional model and the CASE model are in good agreement estimating blow-by. 31 Figure 2.8 Mass Flow Rate at the Outflow Boundary Figure 2.9 Oil Mass Fraction at the Outflow Boundary 32 2.4.4 Reverse Blow-by The vector traces at the inflow boundary, at the selected CA 388 in Figure 2.10 demonstrate traces in the upward direction, which indicates reverse blow-by. During the late expansion and early exhaust stroke, the in-cylinder pressure decreases from its peak value. Due to the pressure difference between the 1st and 2nd land, an upward jet of flow forms at the 1st ring end gap and drives the flow to the combustion chamber. This reverse blow-by or blowback is the amount of gas that blows back from the inter-ring pack to the combustion chamber. It is found to play an important role in determining oil consumption in internal combustion engines [64][65][66]. Figure 2.11 shows the mass flow rate at the inlet boundary throughout the cycle. Positive values indicate flow out of the boundary or backflow, which is assumed to be the reverse blow-by; negative values indicate flow into the boundary or forward flow. Following the same procedure described in the previous section, positive discharged mass per cycle is calculated to be 2.3345e- 05 kg and the mass flow rate out of the boundary or backflow rate is 3.8908e-04 kg/s. CASE estimates the averaged backward flow rate past top ring to be 1.52 e -05 kg/s. Note that this estimate by CASE is a combined flow from the top, gap, and bottom of the top ring which does not capture the three-dimensional ring gap effect, hence the slight difference. 33 Figure 2.10 Vector trace at the Inflow Boundary Showing Reverse Blow-by Figure 2.11 Mass Flow Rate at the Inflow Boundary 34 2.4.5 Pressure Distribution To compare with the pressure distribution found from CASE, the pressure across the length of the piston at the end of the cycle is averaged circumferentially over a cut surface. The pressures at inflow, top land, top groove, second land, second groove, third land, oil groove, and piston skirt are plotted in Figure 2.12 with a piston layout to show the location along with the piston. It is observed that they follow the same trend but the pressure values at every location do not match exactly to each other at 720 CA. This could be attributed to the absence of circumferential effect at the ring end gap in the process of averaging. Figure 2.12 Pressure Across the Piston Length: Comparison Between CASE and Three-Dimensional Simulation at 720 CA 2.4.6 Oil Transport In this section, the results of the oil distribution in the piston ring pack and cylinder liner at selected crank angles across the computational domain are introduced and discussed. Potential driving 35 forces for oil transports between the different parts of the piston ring pack are pressure gradients, shear stress from gas flows, ring/liner, and ring/groove relative motion and inertia forces. The oil moves axially, being driven by the inertia forces resulting from the alternating piston or ring motion. Oil is also observed to move in the circumferential direction primarily by the dragging effect of the blow-by gases. Ring motion inside the groove and the resultant squeezing and pumping effect might trigger both axial and circumferential transport of oil [63]. Note that, inertia force is not modeled in CONVERGE, but the boundary conditions used are obtained from CASE which accounts for it. The initial oil mass fraction in the computational domain is displayed in Figure 2.13. Figure 2.13(a), 2.13(b), 2.13(c) and 2.13(d) show the liner, lands, rings, and grooves, respectively. The cylinder liner is initialized to have oil fraction 0.05, 0.1,0.15 and 1 marked by LR1, LR2, LR3 and Lskirt region of the liner (refer to Table 4 and Figure 3 for initialization). Similarly, the lands, rings, and grooves are also initialized as oil mass fraction to be 0.05 for the top ring, 0.15 for the second ring and 0.15 for the oil ring. The bottom part of the cylinder liner, named Lskirt, and the piston skirt are not displayed for better visualization of the upper part of the cylinder kit consisting of the three rings and adjacent areas. During intake stroke from 0°to 180°CA, the piston moves downward, and scraping occurs during a downward stroke where the oil ring is the leading edge. Figure 2.14 shows the lands, grooves, and rings at the end of the intake stroke in images (a), (c), (d) and (e) respectively. Due to the down-scraping of the taper-shaped second ring, the oil is transported on the 3rd land, which is located below the 2nd ring around BDC (Bottom Dead Center). Figure 14(a) shows that on the third land there are oil puddles more than 0.15 mass fraction. These are colored orange and red. The oil groove shows mass fraction of about 0.18, which is 20% higher than the initialized value. 36 Possibly this oil was driven by the pressure gradient created by the inertia force in the vicinity of the ring/groove clearance transported from the third land. This movement is confirmed by the pressure gradient between the third land and oil groove at 180 CA, as illustrated in Figure 2.14(b). The image shows that the pressure at the third land is 118kPa, whereas in the groove the pressure is 100kPa. The second land does not show any oil mass fraction increase as an effect of top ring scraping. Rather, in the second land, the mass fraction is mostly 0.0675 which is 32.5% less than the initialized value. This indicates that oil must have been transported from the second land to the second groove. In the second groove, as can be seen in Figure 2.14(e), there are oil puddles with mass fraction 0.1125 to 0.1575, which is a 12.5% - 57.5% rise from the initial condition. This means a contribution of one or both of the following phenomena: oil from top ring scraping; and oil from second land transportation to the second groove. Again, Figure 2.14(d) shows that ring scraping creates an oil accumulation on the front side of both the second and the top ring. The top ring has mostly oil mass fraction, 𝑌𝑜𝑖𝑙 = 0.0675 𝑡𝑜 𝑌𝑜𝑖𝑙 = 0.1125; which is a 0.35% to 1.25% rise. The second ring loses some of the initial oil, leading to 𝑌𝑜𝑖𝑙 = 0.0675 in most parts. In some places, it gains oil and rises to 𝑌𝑜𝑖𝑙 = 0.1125. It is observed from Figure 2.14(d) that at the end of the intake stroke, all three rings have 𝑌𝑜𝑖𝑙 higher than the initialized value. At the end of the intake, most of the oil is scraped off from the liner. Then there is only 0.0675 fraction of initialized oil or 32.5% left on LR2. There is oil mass fraction of 0.0225 on the LR1 region, as demonstrated in Figure 2.14(c). In most parts of LR3, the oil mass fraction reduces to 0.0675-0.135, which is a 55%-10% reduction from the initial mass fraction. This also indicates some scraping by the oil ring. Since the LR3 region is adjacent to the piston skirt region where oil is initialized as 100%, some oil flows to the LR3 region; making the oil fraction to be 0.18 in 37 places which are 20% higher than the initialized oil. Figure 2.13 Initial Oil Mass Fraction in the Computational Domain at 0 CA (a) liner (b) lands (c) rings (d) grooves Figure 2.14 Comparative Scenario in the Computational Domain at 180 CA (a) lands (b) pressure gradient between third land and oil groove (c) liner (d) rings (e) grooves 38 Figure 2.15 Oil Mass Fraction across the Domain at 278 CA (a) lands (b) piston skirt (c) pressure gradient through lands (d) rings (e) grooves (f) liner During the compression and early expansion strokes, gases flow from the combustion chamber through the top ring gap into the second land due to the pressure gradient across the top ring. As a result, the pressure in the second land clearance increases because of its finite volume, and a pressure gradient evolves between the second and third land regions. This pressure gradient induces a gas flow from the second land through the second ring gap into the third land. This gas flows most likely along the circumferential direction of the piston's lands. The blow-by gas flow moves the oil on the top land through the ring gap to the second land due to interfacial shear stresses [67]. If an oil control ring with a single gap is used, such as the two-piece oil control ring used in the current study, then the blow-by gas flow on the third land is believed to be mainly in the circumferential direction from the second ring gap to the third ring gap [63]. Figure 2.15 shows a similar trend at 278 CA. Figure 2.15(a) shows that the top land has barely any oil; there is a small oil puddle of mass fraction between 0.025-0.045 on the second land. There 39 is still a small amount of oil in the third land, which indicates that the net oil flow is toward the crankcase. The net effect of inertia force and the dragging action of blow-by gases force most of the oil present in the control volume towards the crankcase region, which is evident in the piston skirt image, Figure 2.15(b), at 278 CA. Figure 2.15(c) shows the gradual pressure gradients across the lands being 400 kPa, 300 kPa, and 200 kPa respectively. Figure 2.15(d) shows the oil mass fraction on the rings at 278 CA. A small puddle of oil of mass fraction Yoil = 0.045 and oil on the second ring imply up scraping. The almost dry LR1 and a small amount of oil on the LR2 region shown in Figure 2.15(f) confirm that. Figure 2.15(e) shows the oil mass fraction in the grooves at 278 CA. The top groove has a small fraction of oil with Yoil = 0.045. A higher amount of oil is accumulated on the second groove. This oil accumulation could be due to one or more of several factors: the combination of inertia, pressure gradient and resulting dragging of gas flow from land to groove; shear-driven oil flow due to the lateral motion of the piston relative to the rings; or pumping due to ring motion [63][14]. The high oil mass fraction on the oil ring groove indicates that most of the oil present is being pushed back toward the crankcase. This movement should persist through the early expansion stroke; but due to lack of continuous supply of oil, by 360 CA there is no oil on the lands, liner regions LR1, LR2, LR3, or the top two rings as displayed in Figure 2.16. The third groove is still immersed in oil with mostly Yoil = 0.18 on the bottom flank of the groove. On the top flank, there are oil puddles of Yoil = 0.0225 to Yoil = 0.135 [see Figure 2.16(c)]. With most of the oil scraped off, the oil fraction keeps dropping gradually throughout the liner. By 360 CA there is almost no oil on LR1 and LR2. There are small collections of oil in places on LR3 because of its adjacency with the skirt region, as can be seen in Figure 2.16(d). 40 During the late expansion and early exhaust strokes, when the gas pressure on the second land (below the top ring) exceeds the pressure in the combustion chamber, reverse blow-by occurs. Some oil from the crankcase or the piston skirt region can be transported from the third land through the second ring end gap, and eventually from the second land through the top ring end gap in the direction of the combustion chamber with the reverse blow-by gas flow. This is depicted in Figure 2.17 and Figure 2.18. Figure 2.16 Oil Mass Fraction across the Piston and Liner at 360 CA (a) lands (b) rings (c) grooves (d) liner 41 Figure 2.17 Oil Mass Fraction across the Piston and Liner at 460 CA (a) lands (b) pressure gradient through lands (c) rings (d) grooves 42 Figure 2.18 Oil Mass Fraction across the Domain at 550 CA (a) lands (b) pressure gradient through lands (c) rings (d) grooves It is observed from Figure 2.17(a) that oil is moving from the third land to the second land through the ring gap. The pressures at the third land, second land, and top land at 460 CA are 600kPa, 500kPa, and 400kPa respectively, as shown in Figure 2.17(b), whereas the in-cylinder pressure is 347.5 kPa. The pressure gradient between the second land and top land is not high enough to drag oil with reverse blow-by gas at that point. At 460CA the top two grooves are almost dry whereas the third groove shows the oil abundance, as illustrated in Figure 2.17(d). At 550 CA, a distinct pressure gradient is observed from Figure 2.18(b). The pressure at third land, second land, and top land are 260kPa, 180kPa, and 100 kPa, respectively. At 550 CA the in- cylinder pressure is 135.360 kPa. This triggers oil to flow from the third land to the second land 43 and then eventually top land through the gap, as shown in Figure 2.18(a). Oil throw-off from the top land has been observed in [12,68]. Oil flow through the top ring gap was observed on a one- cylinder research engine at low load conditions when the top ring was pinned [69,70]. There is some oil in the top two rings which indicates up scraping, and eventual oil accumulation on the other side of the grooves is noticed in Figures 2.18(c) and 2.18(d). At the end of the cycle, some oil is observed in the second and third land, leaving the top land completely dry, as seen in Figure 2.19(a). Figure 2.19(c) shows that the top and second grooves have oil puddles of a very small fraction, with top groove 𝑌𝑜𝑖𝑙 = 0.045 and in the second groove, 𝑌𝑜𝑖𝑙 = 0.045 to 𝑌𝑜𝑖𝑙 = 0.0675. To summarize, the variation of average oil mass fraction in all four regions along the entire piston ring pack is plotted in Figure 2.20(b), with a simplified 2D sketch of the computational domain with different regions identified in Figure 2.20(a). The four images in Figure 2.20(b) depict the following: the top left shows the top ring region, top right shows the 2nd ring region, bottom left shows the oil ring region, and the bottom right one shows the skirt region. A common trend observed in all four images is that there is a drop of oil mass fraction mid-way through the cycle due to the scraping effect during the downstroke. The top ring region dries up first, then the second ring region follows. Oil mass fraction in the oil ring region and piston skirt region does not become completely zero, although a drop is observed. Consequently, the oil mass fraction at the end of the cycle is not the same as it was at the beginning. While a fraction of oil is lost to the crankcase and combustion chamber, the major reason for this decrease lies in the lack of continuous supply of oil, as occurs in actual engine operation. Understanding how oil is supplied to the cylinder kit is crucial to the prediction of oil and gas distribution in different regions of the computational domain. 44 Figure 2.19 Oil Mass Fraction across the Domain at 720 CA (a) lands (b) rings (c) grooves (d) liner Figure 2.20 (a) A 2D sketch of the geometry showing the regions of the computational domain (b) oil fractions in different ring regions 45 2.4.7 Oil Evaporation Oil evaporation from the piston-ring-liner system is also believed to contribute significantly to total oil consumption. Several experimental results indicated that oil evaporation from the liner and piston might contribute substantially to oil consumption [56,71–73]. Besides, several purely theoretical approaches studied oil evaporation from the liner and found sensitivities in the evaporation process to oil composition and cylinder liner temperatures [66,74,75]. Qin et al. [76] found oil evaporation to be extremely sensitive to the liner temperature, positively correlated with engine speed, moderately related to liner material but insensitive to oil film thickness. The evaporation of oil from different regions of the engine, such as the cylinder liner, piston, and the oil sump, contributes to total oil consumption. The dominant contribution is believed to be the evaporation from the piston-ring-liner system, as the oil present in these regions is exposed to higher gas flow rates and temperatures [13]. Figure 2.21 Oil Evaporation in the Lands throughout the Cycle (a) 180 CA (b) 270 CA (c) 360 CA (d) 460 CA (e) 550 CA (f) 720 CA (g) A 2D sketch of the geometry showing the regions of the computational domain 46 Figure 2.22 Oil Evaporation in the LR1, LR2 and LR3 Boundaries of the Liner throughout the Cycle (a) 180 CA (b) 270 CA (c) 360 CA (d) 460 CA (e) 550 CA (f) 720 CA (g) A 2D sketch of the geometry showing the regions of the computational domain Figure 2.23 (a) A 2D sketch of the geometry showing the regions of the computational domain (b) oil vapor fractions in different ring regions 47 Oil evaporation across the piston lands throughout the cycle is illustrated in Figures 2.21(a) through Figure 2.21(f), and the regions are identified in a two-dimensional sketch in Figure 2.21(g). It is observed from Figure 2.21(a) that there is a small amount of oil vapor on the top two lands at the end of the intake stroke. After that, oil evaporation is not prominent during the compression stroke; refer to Figure 2.21(b) and to Figure 2.21(c), which indicates the negligible significance of oil evaporation from the piston ring pack into the blow-by flow. During late expansion and early exhaust stroke, oil is transported as vapor into the second land and subsequently to the top land as a result of reverse blow-by, as shown in Figure 2.21(d) and Figure 2.21(e). The oil vapor in the top land could be transported to the combustion chamber. Figure 2.22(a) through Figure 2.22(f) demonstrates oil evaporation in the LR1, LR2, LR3 boundaries of the liner throughout the cycle, and the regions are identified in a two-dimensional sketch in Figure 2.22(g). The Lskirt is not displayed for better visualization of the other three boundaries. It is observed that the likelihood of evaporation on the liner is directly related to the presence of oil. Throughout the cycle, most of the evaporation is observed in the LR3 region because of the high mass fraction of oil present there. There is a very small amount of evaporation noticed on the LR1 and LR2 region during the first two strokes, as can be seen in Figures 2.22(a), 2.22(b) and 2.22(c). During late expansion and early exhaust stroke, oil vapor is observed on both LR1 and LR2 regions; refer to Figure 2.22(d) and Figure 2.22(e). The mass fraction of oil evaporation throughout the cycle is in the minuscule range: on the order of 10−11. Figure 2.23 shows the oil vapor fraction in the four regions. In the present study, 𝐶7 𝐻16 is assumed to be the most volatile constituent of the oil, and the figures above capture the evaporation of 𝐶7 𝐻16 .Because of the adequate amount of oil, the oil vapor fraction is on the order of 10−6 in the oil ring region and skirt region. In contrast, oil seldom evaporates in the top and second ring 48 regions. 2.5 Conclusions The current study is the first modeling effort toward developing a three-dimensional, multi-phase CFD methodology to investigate liquid oil and gas transport in the cylinder-kit assembly. A 1D model is developed, using CASE to generate the boundary conditions of the three-dimensional model. A three-dimensional, multi-phase, CFD-model is developed using CONVERGE. The key findings of the model show oil and gas mass distribution in all parts of the domain and their impact on blow-by, reverse blow-by and oil consumption. This model will be further extended to explore the effects of the various ring design parameters on blow-by and oil consumption. A mesh sensitivity study indicated the need for a very fine mesh size (less than 100 µm) for capturing the oil-air interface as well as certain areas where the ring-groove or ring-liner clearance is considerably low. A study with such finer mesh size would require more computational resources and time. The circumferential flow across the lands is displayed in three dimensions, with consideration of the ring end gap. The ring end gaps and their relative locations are observed to be a key factor behind flow variation across the regions. Comparison of blow-by and pressure distribution along piston length in the quasi-one-dimensional model and the three-dimensional model showed OK agreement. The flow at the inflow boundary gives an estimate of reverse blow-by, which is also compared with CASE results. The results show the distribution of oil in terms of mass fraction during an engine cycle. Mechanisms behind the oil flow patterns are identified qualitatively. Insights from the three- 49 dimensional model can lead to better calibration of the one-dimensional CASE model. Evaporation in the liner and top land is characterized. It is observed that the possibility of liquid oil changing to oil vapor is higher in the regions where oil is distributed in a higher amount. It is observed that continuous oil supply management significantly dictates the estimate of oil lost to the crankcase and combustion chamber. To this end, the oil will be re-initialized in the piston skirt region just before the end of every upstroke in future works. The future works will also incorporate the ring dynamics and its effect on liquid oil-gas transport. 50 Chapter 3 THE EFFECT OF RING-GROOVE GEOMETRY ON ENGINE CYLINDER-KIT ASSEMBLY USING THREE-DIMENSIONAL MULTIPHASE PHYSICS-BASED MODELING METHODOLOGY-PART II 3.1 Abstract Cylinder-kit tribology has been a significant focus in developing internal combustion engines of lower emission, reduced friction and oil consumption, and higher efficiency. This work addresses the impact of ring-groove geometry on oil (liquid oil and oil vapor) transport and combustion gas flow in the cylinder kit, using a dynamic three-dimensional multiphase modeling methodology during the four-stroke cycle of a piston engine. The ring and groove geometry, along with the temperature and pressure conditions at the interface between piston and liner, trigger the oil and gas (combustion gases and oil vapor) transport. A study of the second ring dynamics is presented to investigate the effect of negative ring twist on the three-dimensional fluid flow physics. The oil (liquid oil and oil vapor) transport and combustion gas flow processes through the piston ring pack for the twisted and untwisted geometry configurations are compared. Note that the twisted second ring induces a change in crevice volume as well as the leakage area. This, in turn, affects the overall mass flow rate and influences the liquid oil-gas (combustion gases and oil vapor) distribution in different regions of the cylinder kit. A comparison with the untwisted geometry for this cylinder- kit shows that the twisted second ring resulted in a higher blowby but lower reverse blowby and oil consumption. The oil mass fraction distribution pattern across the cycle also confirms the fact. The comparison of the model predicted oil consumption with existing literature shows that oil consumption is within the reasonable range for typical engines. Implementing this three- dimensional methodology leads to a better understanding of cylinder-kit fluid flow physics and 51 the effect of ring-groove geometry on it. The findings presented in this work pave the way to further the ongoing development of optimum cylinder-kit designs with controlled gas leakage, low oil consumption, and low cylinder-kit friction. 3.2 Introduction The performance of a modern engine cylinder-kit assembly depends on the tradeoff between two major functions- sealing and lubrication. The proper balance of these can ensure durability, emission reduction, and fuel economy improvement. Thus, a great deal of attention is focused on improving the piston ring pack design to control blowby, friction, wear, and oil consumption to address these issues. The ring dynamics, i.e., radial, and axial motion and twist, have significant impact on ring stability, sealing capability, inter-ring pressure, and overall ring performance. Most rings are cut with a bevel to help the ring twist in its groove. Due to the asymmetrical cross- section, when the ring is fitted in the cylinder bore, it is subjected to the three-dimensional distortions called twist. If the ring has a bevel on the inside top edge, the twist is considered to be the positive twist. It causes the ring face to twist in the upward direction toward the piston crown. It allows the cylinder pressure to get behind the ring and push it out against the cylinder wall for a better seal. Positive twist also holds the inside bottom edge of the ring against the ring groove for oil control purposes; see Figure 3.1(a). If the ring has a bevel on its inside bottom edge, it causes the ring face to twist in the downward direction toward the piston skirt. This is considered to be the negative twist. This helps position the bottom face of the tapered edge of the ring against the cylinder liner, while the inside top edge seals against the ring groove. That allows the oil scraped off the cylinder wall to pass under the ring and through the return holes in the back of the ring 52 groove; see Figure 3.1(b) for a negatively twisted tapered ring. The twist angle of a ring is not uniform across the circumference of the ring [77,78]. Cheng et al. [79] showed that for a positive static twisted second ring, the twist angle varies from a small value at the ring end-tip to a larger value at the back point (opposite to the ring end-gap). Figure 3.1 Twist on ring (a) Positive twist (b) Negative twist Ruddy et al. [6,80] addressed the importance of the contact point between rings and their grooves. They found some interesting phenomena related to ring twist, such as possible ring flutter if the ring's contact point and its groove are at the outside edge. Keribar et al. [81] developed an integrated ring pack model considering the effect of the forces between the rings and their grooves. The flow through the ring grooves was, however, simplified as an orifice flow. Tian et al. [82] studied the effect of ring twist and developed a ring dynamics and gas flow model to study ring/groove contact characteristics, ring stability, and blowby. They found ring flutter to occur for the second ring with a negative static twist under normal operating conditions. Ejakov et al. [35] 53 presented a three-dimensional (3D) beam model of a piston ring to explore the influence of 3D ring twist on ring performance and blowby. Mittler et al. [77] developed a 3D analytical tool to analyze ring twist and ring movement under dynamic loading conditions. Cheng et al. [39,79] developed a 3D finite element analysis (FEA) ring contact model integrating ring-cylinder bore- groove interaction, in which they studied the effect of ring twist angle on ring dynamics and gas dynamics. Of all the published works, the effects of ring twist on cylinder-kit tribology demonstrating three-dimensional multiphase (liquid oil, combustion gas, oil vapor) fluid transport across the cycle have not been documented in 3D before. The present study utilizes the three-dimensional multiphase physics-based modeling methodology presented by Chowdhury et al. in [83]. In Part I [83], the details of the model development methodology have been outlined to explore liquid oil-gas (combustion gases and oil vapor) transport inside the cylinder-kit assembly. The work here presents the effect of ring twist on the overall liquid oil-gas (combustion gases and oil vapor) transport along the piston ring pack in a three-dimensional model for the first time. A negative static twist value is applied to the second compression ring, and the effect of dynamically varying twist geometries over the cycle is studied. The ring dynamic twist can introduce circumferential variations to ring/liner lubrication and influence the overall performance of the ring-pack in terms of friction and oil consumption. 54 3.3 Workflow Figure 3.2 From tools to the solution (workflow) The entire workflow of the model development consists of four distinct stages outlined as Block 1, Block 2, Block 3, and Block 4 in Figure 3.2.In the first stage, following the procedure as described in Part I [83], a quasi-one-dimensional model is developed in CASE [29], applying a negative static twist of -0.15 to the taper-shaped second compression ring. Primary engine dimensions and operating conditions are presented in Table 3.1. Cylinder-kit geometrical dimension, material properties, flow characteristics, thermodynamic attributes are some key input parameters for the model development in CASE, as depicted in Block 1 of Figure 3.2. Inter-ring pressure and crank angle resolved piston-ring motion from the CASE analysis is exported to be the boundary conditions for the second stage, i.e., 3D ring FEA contact model. 55 Table 3.1 Engine operating conditions [83] Parameter Value Bore 50.6mm Stroke 50mm Connecting rod length 130mm Engine speed 2000 rpm Compression ratio 14 Engine oil SAE 10W-30 In the second stage, as shown in Block 2 of Figure 3.2, the ring cross-section geometry, ring material properties, ring outer diameter (OD) profile at its free state (known as ring free shape profile); inter-ring pressure, ring velocity, and ring acceleration obtained from CASE are fed into a three-dimensional ring FEA contact model developed by Cheng et al. [39,79]. Note that for the second ring, inter-ring pressure means the pressure at the adjacent boundaries, namely- 1. Second land- This is the surface between the top and second groove 2. Second groove- The taper faced second compression ring is inserted in this groove 3. Third land- This is the surface between the second and third groove 56 Figure 3.3 Inter-ring pressure for second ring input to 3D FEA contact model generated by CASE[29] The pressure across the cycle at these three boundaries with a schematic showing the boundaries is presented in Figure 3.3. The in-cylinder pressure is superimposed for reference. The sharp peaks between 290 CA and 442 CA in Figure 3.3 indicate that the ring flutters. The pressure across the cycle at these three boundaries with a schematic showing the boundaries is presented in Figure 3.3. The in-cylinder pressure is superimposed for reference. The sharp peaks between 290 CA and 442 CA in Figure 3.3 indicate that the ring flutters. This means an unstable axial in-groove motion, which is the result of strong coupling between the ring dynamics (ring inertial force) and flow dynamics (pressure force). Ring flutter can change the stability of the second ring and flow paths. Tian et al. [82] observed second ring flutter due to a negative static 57 twist during the late compression and early expansion strokes [82]. Cheng et al. [84] studied the second ring dynamics and concluded that the second ring is more likely to flutter in the presence of a negative static twist. The second ring flutter phenomenon was also observed from the experiment by Furuhama et al. [85] for various operating conditions. Figure 3.4 Second ring motion across the cycle input to 3D FEA contact model (a) ring velocity (b) ring acceleration The other inputs to the contact model are ring velocity and acceleration, demonstrated in Figure 3.4(a) and Figure 3.4(b), respectively. The ring dynamics over an engine cycle is analyzed using the finite element method with eight-node hexahedral elements by employing the penalty method [86]. The 3D analytical ring model [39,79] addresses the pressure/force distribution and ring conformability between the ring-cylinder bore interface and ring-groove interface. The contact pairs between ring-liner and ring-groove interface are computed via the penalty method [86] and 58 the force release technique. For each crank angle, the twist angle variation over circumferential locations, deformed ring shape profile, piston ring OD profile, and ring location is calculated. The free shape ring mesh and the deformed ring mesh generated using the contact model [39,79] are shown in Figure 3.5. Figure 3.5(a) shows the three-dimensional view. A two-dimensional image is presented in Figure 3.5(b) for better visualization. The cyan mesh in Figure 3.5 represents the free shape ring, while the purple mesh represents the deformed ring shape. It is evident that the ring is pushed inward from its free state. It is hard to differentiate the change in geometry with the crank angle from the ring mesh of Figure 3.5. For this purpose, the twist angle vs. circumferential location plot at each crank angle is carefully inspected to find out how the ring twist angles are affecting the ring geometry dynamically. Figures 3.6(a), 3.6(b), 3.6(c), and 3.6(d) show the general trend of the twist angle variation at each stroke, respectively, for one half of the ring. The other half is assumed to be symmetric. The ring twists negatively from the ring back to the ring end. From the ring back, the twist angle first reduces and then increases toward the ring end. Five circumferential locations are picked to characterize the twist angle variation, namely 1, 45, 88, 131, and 168 degrees. These are denoted by A, B, C, D, and E, respectively, in Figures 3.6(a), 3.6(b), 3.6(c), and 3.6(d). Note that 0 degree refers to the ring back, and 168 degree refers to the ring end gap. The circumferential locations are shown in a schematic in Figure 3.6. After close observation of the variation of the twist angle over the circumference, twenty-four (24) crank angles are selected. The twist angle variation is significant at those 24 crank angles. This means the ring location and ring outer diameter (OD) profile undergo a transformation as well throughout the cycle. The selected 24 crank angles at which the ring location and OD change 59 significantly are listed in Table 3.2. Because of ring flutter between 290 CA and 442 CA (see Figure 3.3), the ring location and the ring OD profile change frequently. This, in turn, results in a frequent change in the cylinder-kit geometry during this crank angle duration. It is observed from Table 3.2 that nineteen (19) of the total 24 crank angles correspond to that late compression-early expansion crank angle interval. So, it is evident that the geometry changes frequently due to ring flutter induced by the negatively twisted ring. Figure 3.5 Free shape and deformed ring mesh generated using the ring contact model [39,79] (a) three- dimensional view (b) two-dimensional view 60 Figure 3.6 General trend of twist angle vs. circumferential location: (a) Intake stroke from 0 CA to 180 CA (b) Compression stroke from 181 CA to 360 CA (c) Expansion stroke from 361 CA to 540 CA (d) Exhaust stroke from 541 CA to 720 CA The ring location and OD profile at the above-mentioned 24 crank angles are then fed to the next stage LINCC [40], as shown in Block 3 of Figure 3.2. As mentioned earlier, the geometry changes 24 times throughout the cycle. Those twenty-four dynamically changing geometries are generated using the ring location and OD profiles from the contact model and other geometrical attributes of the cylinder kit. Figure 3.7 shows a section cut view of one such geometry in Stereolithographic (STL) format. The yellow arrows indicate the generic flow path. The flow domain in the created geometries would be solved in the next stage. Therefore, these dynamically changing geometries from LINCC are used to develop a three-dimensional multiphase model in CONVERGE [46], as shown in Block 4 of Figure 3.2. 61 Table 3.2 Selected crank angles at which the ring location and OD change CA CA CA CA 0 331 372 406 181 353 377 415 276 361 382 420 310 365 387 432 321 367 398 458 326 369 401 541 Figure 3.7 Section cut view of the geometry generated using LINCC [40] The simulation starts at 0 CA with the first geometry and does not change until 181 CA. At 180.99 CA, the results are mapped as the initial condition for the second geometry, and the process goes on throughout the cycle. 62 The quality of the generated mesh is an important parameter in solving the flow physics of Internal Combustion (IC) Engines. The previous research studies show that the mesh structure and grid accuracy significantly affect the quality of flow solvers, especially for complex configurations. (see e.g., [87–89]). The grid sensitivity study that we conducted in [83] showed that the grids need to be sufficiently fine to resolve the oil-air interface and to capture the tiny ring-groove clearances on the micron scale. Based on that study, the grid size is set at 75 µm in x and y directions and 100 µm in the z-direction, resulting in 4,196,137 cells in the computational domain, see Figure 3.8. The mesh is generated automatically in CONVERGE using a cut-cell cartesian meshing technique resulting in a body-fitted volume mesh. The interior cells are hexahedrals. The cells intersected by the geometry surface are cut by the surface to form arbitrary-sided polyhedra [46].It required 760,000 CPU hours to complete the simulation. Figure 3.8 Generated mesh of the domain in CONVERGE [46] To analyze the liquid oil-gas multiphase flow, the volume of fluid model (VOF) [48] is used. It is assumed that the gas phase is a combination of N2 and O2 while SAE 10W30 oil is the liquid phase. The homogeneous relaxation model [57] is used to model the evaporation, assuming that 63 n-heptane (C7 H16 ) is the most volatile constituent of the liquid oil. For liquid oil- combustion gas mass fraction initialization, the domain is divided into four regions; more details can be found in a later section of the paper. The in-cylinder pressure and temperature boundary conditions across the domain are shown in Figure 3.9 and Figure 3.10, respectively; see [83] for the detailed description of the methodology including boundary conditions and governing equations. Figure 3.9 Schematic View of the Boundary and Initial Conditions [83] 64 Figure 3.10 In-cylinder Pressure Profile of the Engine at 2000 rpm [83] 3.4 Results & Discussion 3.4.1 Circumferential flow To compare the general circumferential flow behavior, the stream traces on the three lands of the untwisted and the twisted cylinder-kit geometry at a selected CA (Crank Angle) 369 is depicted in Figure 3.11, Figure 3.12, and Figure 3.13, respectively. It is observed that the flow velocity across the lands is influenced by the ring end gap in both untwisted and twisted configuration. In the top land, flow is mostly axial, but small differences are observed, near the right-hand side, between the untwisted configuration [Figure 3.11(a)] and the twisted one [Figure 3.11(b)]. Figure 3.12(a) and Figure 3.12(b) shows the flow in the second land with a view of the top ring end gap and second ring end gap for the untwisted and twisted configuration, respectively. The flow in the second land is greatly affected by the respective location of the first and second ring end gaps. 65 Small vortices ranging between 13 m/s-25 m/s are formed just below the top ring end gap in the twisted geometry, see Figure 3.12(b), whereas the stream traces range between 5 m/s-19 m/s at the same location of the untwisted geometry [Figure 3.12(a)]. Above the second ring end gap, the flow velocity is mostly 5 m/s-11 m/s for both configurations with small difference, as shown in Figure 3.12(a) and Figure 3.12(b). It is clear that the differences observed could cause local temperatures and oil evaporation rates to change between the twisted and untwisted configurations. Figure 3.11 Top land Stream traces at a selected crank angle (369 CA) (a) Untwisted configuration (b) Twisted configuration 66 Figure 3.12 second land Stream traces at a selected crank angle (369 CA) (a) Untwisted configuration (b) Twisted configuration. Figure 3.13 Third land Stream traces at a selected crank angle (369 CA) (a) Untwisted configuration (b) Twisted configuration. 67 The flow circulating structures exiting the second ring end gap into the third land are larger than on the first or second lands. One vortex formed on each side of the 2nd ring end gap, but the gradients in the complicated structures differ between the untwisted [Figure 3.13(a)] and the twisted [Figure 3.13(b)] configurations. The flow velocity reaches nearly 120 m/s below the second ring gap for both configurations. More stream traces with higher velocities are seen in the untwisted configuration [Figure 3.13(a)] compared to the twisted one [Figure 3.13(b)]. Again, these different flow fields will result in variations in heat transfer and oil consumption, as shown in the following sections. 3.4.2 Estimating Blow-by To calculate the blowby, mass flowrate at the outflow boundary is plotted in Figure 3.14. In this study, the engine runs at 2000 RPM, and thus it takes 0.06 seconds to complete a cycle. Therefore, 1 1 CA corresponds to 12 millisecond. The discharged mass in kg is found by integrating the area under the mass flow rate versus time graph. The discharged mass is calculated to be 1.8441e-04 kg. This cumulative mass is the sum of the mass of oil and gas (combustion gases and oil vapor). Observing the oil mass fraction at the outflow boundary in Figure 3.15, the averaged liquid oil, combustion gas, and oil vapor mass fraction are found to be 0.88, 0.115, and 1.56e-05, respectively, at the outlet boundary. Therefore, the discharged combustion gas per cycle is found to be 3.56e-04 kg/s. The final step to calculate blowby is to convert the mass flow rate from kg/s to liter/min. It is assumed that air density is 1.1766 kg/m3 considering air pressure and temperature to be 101325 Pa and 300K, respectively. Finally, the blowby with the twisted second ring is calculated to be 18.4 liters/min. The untwisted geometry results in a blowby of 9.05 liters/min. So, it appears that 68 the second ring twist causes a rise in the blowby at this operating condition. This can be attributed to the flutter caused by the pressure and inertia loading on the negatively twisted second ring. This means a higher amount of combustion gas flow is released toward the crankcase, resulting in higher blowby. Figure 3.14 Mass flow rate at the outflow boundary for twisted configuration 69 Figure 3.15 Oil mass fraction at the outflow boundary for twisted configuration 3.4.3 Estimating Reverse Blow-by & Oil Consumption from the Cylinder Kit Estimation of reverse blowby and oil consumption requires the measurement of the amount of combustion gas (from the inter-ring pack) and liquid oil flowing (from the crankcase through the ring pack) out of the inflow boundary toward the combustion chamber during a complete cycle. These are calculated by summing the mass flow rate that leaves the inflow boundary over the whole cycle. Figure 3.16 shows the mass flow rate at the inflow boundary throughout the cycle. Positive values indicate flow out of the boundary or backflow, which is assumed to be the reverse blowby; negative 70 values indicate flow into the boundary or forward flow directed toward the crankcase. Following the same procedure described in the previous section, the positive discharged mass per cycle (backflow or reverse flow) is calculated to be 2.7859e-05 kg, and the negative discharged mass per cycle (forward flow) is calculated to be 4.25840e-05 kg. This means that the discharged mass toward the combustion chamber is 34.57% lower than that toward the crankcase via the inflow boundary. Since at the inflow boundary, we are concerned about the reverse blowby and oil consumption calculation, the positive discharged mass of 2.7859e-05 kg is taken for further investigation. This cumulative mass is the sum of the mass of oil and mass of gas (combustion gases and oil vapor). Observing the oil mass fraction at the inflow boundary in Figure 3.17, the averaged combustion gas, liquid oil, and oil vapor mass fraction are found to be 0.99, 0.00102, and 1.22e-12, respectively. Figure 3.16 Mass flow rate at the inflow boundary for twisted configuration 71 Figure 3.17 Oil mass fraction at the inflow boundary for twisted configuration Thus, the mass flow rate out of the boundary or backflow rate is 4.63845e-04 kg/s. This results in a reverse blowby of 23.99 liters/min for the twisted second ring simulation, whereas, for the untwisted geometry, the reverse blowby is found to be 24.34 liters/min. Thus, the twisted second ring reduces the reverse blowby by a small margin. Similarly, the liquid oil flowing toward the combustion chamber per cycle is calculated to be 0.47 mg/s or 1.7 g/hr for the twisted second ring configuration. This is found to be 0.667 mg/s or 2.4 g/hr for the untwisted geometry. Typical engine oil consumption is 2 mg/s per liter of engine displaced volume [90]. For the 0.1 L engine operated at 5.3 bar IMEP and 2000 rpm, the ballpark 72 estimate of oil consumption is 0.2 mg/s or 0.72 g/hr. Thus, the model predicts three (3) times and 2.4 times higher than the typical value for the untwisted and the twisted configurations, respectively. Engine oil consumption can vary up to ten (10) times, depending on the load-speed variation and operating conditions [13]. Kraus et al. [91] reported that oil consumption increases with engine size, and Shore et al. [92] confirmed the variation of oil consumption with engine size. Accordingly, the model prediction is within the range for typical engines. Table 3.3 Selected specifications of the engines used in literature Author No. of Cylinder Oil Used Displacement (L) Froelund et al. V6 SAE 10W30 3.8 Ariga et al. 6 inline SAE 10W30 7.6 Apple et al. 4 inline SAE 5W40 1.9 Delvigne et al. 4 inline SAE 5W30 2 Soejima et al. 4 inline SAE 10W30 2.347 Wong et al. 1 SAE 10W30 0.309 Current work 1 SAE 10W30 0.1 The oil consumption is also compared with experimental oil consumption measurements from existing literature [93–99] where the engine is operating at 2000 RPM in Figure 18. The researchers in [93–99] employed tracer methods to measure oil consumption. Wong et al. [98] used tritium tracer for a 0.3L single-cylinder Kubota EA300N engine. Ariga et al. [94] conducted an oil consumption measurement system using the Sulfur-Di-Oxide (SO2 ) tracer method under various engine operating conditions, including a 200-hour durability test. Froelund et al. [93] measured oil consumption using SO2 analyzers technique, where the engine is operated on a sulfur-containing lubricant and a zero-sulfur fuel. Soejima et al. [97] also used sulfur tracer. Apple et al. [95] presented a high-temperature oxidation method (HTOM) to estimate engine oil 73 consumption in real-time. Atomized lubrication oil and diesel engine exhaust were used to evaluate the HTOM performance. Delvigne et al. [96] carried out a study based on the use of Germanium- 69 (Ge-69), a radioisotope that substitutes carbon atoms in base oil molecules. They studied the impact of lubricant formulation on oil consumption. Several key specifications of studies [93–99] are listed in Table 3.3. Figure 3.18 shows the oil consumption in g/L-hr vs. brake mean effective pressure (BMEP) in bar for the studies mentioned. The red star and the green star denote the oil consumption predicted by the model for untwisted geometry and twisted configuration, respectively. The comparison shows that the oil consumption predicted by our model is considerably higher than those reported in the literature. Again, there are experimental variability, engine to engine or ring pack to ring pack design, and performance variation that might have attributed to this difference to some extent. Another key takeaway of the model prediction is that the twisted second ring has slightly reduced oil consumption. The second ring flutter, induced by the twist, causes a higher blowby. More oil is triggered down to the third land and eventually crankcase. This, in turn, results in lower oil flow being entrained toward the combustion chamber. 74 Figure 3.18 Oil consumption vs. brake mean effective pressure at 2000 rpm; comparison with literature 3.4.4 Comparison of Oil Mass Fraction Distribution Along the Piston Length between twisted and untwisted geometry We introduced and discussed the results of the liquid oil and oil vapor distribution in the untwisted piston ring pack and cylinder liner at selected crank angles across the computational domain in detail in [83]. We described the comparison among land-groove pressures with the in-cylinder pressure and the resultant effect on liquid oil mass fraction distribution and phenomena like blowby or reverse blowby. It is observed that the general trend remains the same across the cycle for the twisted configuration. 75 Figure 3.19 A 2D sketch of the geometry showing the regions of the computational domain. In this section, the comparison of the liquid oil mass fraction between the twisted and untwisted configuration is presented. The pressure gradient among different regions of the 3D model is not shown in detail to avoid redundancy from [83]. 76 Figure 3.20 Initial oil mass fraction on the lands at 0 CA for both twisted and untwisted configuration At the beginning of the cycle, for both twisted and untwisted configurations, liquid oil mass fraction 𝑌𝑜𝑖𝑙 is initialized to be 0.05, 0.1, 0.15, and 1 at the top ring region, second ring region, oil ring region, and skirt region, respectively (details in [83]). Figure 3.20 shows the initial distribution along the three lands of the piston at 0 CA. The lands are shown in a 2D sketch of the geometry in Figure 3.19. At the end of the intake stroke, i.e., 180 CA, the oil transport effect due to ring down scraping and pressure gradient is higher on the untwisted geometry. For the twisted geometry, the third land has 𝑌𝑜𝑖𝑙 ranging between 0.1125 and 0.1575, compared to the untwisted one where 𝑌𝑜𝑖𝑙 is 0.045-0.09. On the second land, the untwisted geometry has 𝑌𝑜𝑖𝑙 = 0.0225, whereas the twisted geometry shows oil puddles of 𝑌𝑜𝑖𝑙 = 0.045. The oil mass fraction on the untwisted and twisted geometry lands are shown in Figure 3.21(a) and Figure 3.21(b), respectively. 77 Figure 3.21 Oil mass fraction on the lands at 180 CA (a) untwisted configuration (b) twisted configuration. Figure 3.22 Oil mass fraction on the lands at 278 CA (a) untwisted configuration (b) twisted configuration. 78 Figure 3.23 Pressure on the lands at 278 CA (a) untwisted configuration (b) twisted configuration Figure 3.24 Oil mass fraction on the lands at 550 CA (a) untwisted configuration (b) twisted configuration 79 Figure 3.25 Pressure on the lands at 550 CA (a) untwisted configuration (b) twisted configuration Figure 3.26 (a) Liquid oil mass fractions in different ring regions comparison between twisted and untwisted configuration (b) A 2D sketch of the geometry showing the regions of the computational domain. During the compression and early expansion strokes, gases flow from the combustion chamber through the top ring gap into the second land due to the pressure gradient across the top ring. As a 80 result, the pressure in the second land clearance increases because of its finite volume, and a pressure gradient evolves between the second and third land regions. This pressure gradient induces a gas flow from the second land through the second ring gap into the third land. This gas flows most likely along the circumferential direction of the piston's lands. The net effect of inertia force and the dragging action of blowby gases force most of the oil present in the control volume towards the crankcase region [83]. It is observed that oil flowing toward the crankcase due to the blowby effect is higher in mass fraction for the twisted geometry. At 278 CA, the oil mass fraction on the third land of the twisted geometry is 𝑌𝑜𝑖𝑙 = 0.1125 and higher, as seen in Figure 3.22(b) compared to Figure 3.22(a). The untwisted geometry in Figure 3.22(a) shows, the oil mass fraction on the third land is between 0.045 and 0.09. At 278 CA, the pressure across the top, second, and third land for the untwisted configuration is 300 kPa, 289 kPa, and 151 kPa, respectively, shown in Figure 3.23(a). For the twisted configuration, the pressure across the lands is 300 kPa, 298 kPa, and 128 kPa, respectively; see Figure 3.23(b). During the late expansion and early exhaust strokes, when the inter-ring pressure exceeds the pressure in the combustion chamber, reverse blowby occurs. Some oil from the crankcase or the piston skirt region can be transported from the third land through the second ring end gap, and eventually from the second land through the top ring end gap in the direction of the combustion chamber with the reverse blowby gas flow [83]. It is observed that the oil entrained by the reverse blowby reduces significantly due to the twisted 2nd ring. Figure 3.24 shows the comparison at 550 CA between the untwisted geometry [Figure 3.24(a)] and the twisted geometry [Figure 3.24(b)]. There are oil puddles on the second land and top land on the untwisted geometry, which are not present on the twisted geometry. For the untwisted configuration, the pressures at the third land, second land, and top land at 550 CA are 260 kPa, 157 kPa, and 135.4 kPa respectively, as shown 81 in Figure 3.25(a). For the twisted configuration, the pressures at the third land, second land, and top land at 550 CA are 290 kPa, 156 kPa, and 137 kPa, respectively, as shown in Figure 3.25(b). The in-cylinder pressure is 135 kPa. To summarize, the variation of average oil mass fraction in all four regions along the entire piston ring pack for both untwisted and twisted geometry is plotted in Figure 3.26(a). A simplified 2D sketch of the computational domain with the regions labeled is shown in Figure 3.26(b). The four images in Figure 3.26(a) depict the following: the top left shows the top ring region, the top right shows the 2nd ring region, the bottom left shows the oil ring region, and the bottom right one shows the skirt region. It is observed that, at the end of the cycle, the oil mass fraction is lower in the top and second region and higher in the oil ring region and skirt region for the twisted geometry. The twisted second ring induces a change in crevice volume as well as leakage area. This, in turn, affects the overall mass flow rate and influences the gas-oil distribution in different regions of the cylinder kit. The higher oil mass fraction in the regions below the second ring indicates the net higher oil flow toward the crankcase. The oil mass fraction pattern across the cycle confirms the increase in blowby and reduction in reverse blowby and oil consumption for the twisted geometry. 3.5 Conclusions The current study investigates the effect of ring twist on cylinder-kit geometry and tribology, utilizing a novel three-dimensional, multi-physics methodology developed by the authors. This is the first attempt to demonstrate the effect of ring twist on the overall liquid oil-gas (combustion gases and oil vapor) transport along the piston ring pack in 3D. A 1D model is developed, using CASE to generate the boundary conditions of the three-dimensional model. A 3D ring FEA contact 82 model generates the circumferential variation of the twist angle and subsequent variation in ring location and ring OD profile. These, together with cylinder-kit cross-section geometry, are used to generate the dynamically varying geometry across the cycle. Finally, a three-dimensional, multiphase, CFD-model is developed using CONVERGE. Note that instead of using the ring motion as a boundary condition for the 3D CFD model, the dynamically changing geometry configurations are superimposed throughout the simulation. The model's key findings show the effect of second ring twist on liquid oil-gas distribution throughout the domain and on properties such as blowby, reverse blowby, and oil consumption. It is observed that the second ring negative twist induces a change in the crevice volume and leakage area. The resultant change in the mass flow rate affects the liquid oil-gas distribution in the domain. A comparison with the untwisted geometry shows that the twisted second ring results in a higher blowby but lower reverse blowby and oil consumption. Note: the results are generated to demonstrate the capabilities of the methodology. They would be different in other geometries, contingent on the specific system geometry and operating condition. The methodology and the findings presented in this work can play a significant role in cylinder-kit performance assessment and design optimization. To this end, we are implementing the methodology on a full-size real engine for experimental validation which will be demonstrated in future. 83 Chapter 4 ONE DIMENSIONAL OPTIMAZATION STUDY OF CUMMINS ACADIA ENGINE 4.1 Introduction Cylinder-Kit Analysis System for Engines (CASE) 1-D model is used both as a tool for existing design assessment and for experimenting with new cylinder-kit design ideas. For users to have confidence in the results of simulation studies, it is essential to make sure that the models are properly evaluated. Calibration is a major element to this evaluation and refers to the estimation and adjustment of model parameters to improve the agreement between model output and experimental results [100]. In CASE, there are ten parameters related to flow paths that need to be configured and adjusted to match with the experimental blowby and land pressure results. These parameters often have nonlinear effects on the model’s behavior which makes it difficult to predict how the model will behave under new parameter configurations. Hence, the process of manipulating a model’s parameters to match experimental conditions is often nontrivial. Once a quantitative measure of similarity between the model data and experimental results are developed, optimization methods [101] are used to explore the parameter space of the model in search of the optimal configuration. In this chapter, an optimization study using CASE [29] and HEEDS [102]. Both programs are linked together for the study where HEEDS reads through the input and output files from CASE to give an output for the best design that meets the requirements of the study. HEEDS [102] uses a hybrid and adaptive algorithm called SHERPA [103] as its default search method. SHERPA is a proprietary hybrid and adaptive search strategy available within the HEEDS software code. SHERPA means “Simultaneous Hybrid Exploration that is Robust, Progressive, and Adaptive”. It employs multiple search strategies at once and adapts to the problem as it “learns” about the design space. It is a direct optimization algorithm in which all 84 function evaluations are performed using the actual model as opposed to an approximate response surface model. During a single parametric optimization study, SHERPA uses the elements of multiple search methods simultaneously (not sequentially) in a unique blended manner. This approach attempts to take advantage of the best attributes of each method. Attributes from a combination of global and local search methods are used, and each participating approach contains internal tuning parameters that are modified automatically during the search according to knowledge gained about the nature of the design space. This evolving knowledge about the design space also determines when and to what extent each approach contributes to the search. In other words, SHERPA efficiently learns about the design space and adapts itself to effectively search all sorts of design spaces, even very complicated ones. All the parameters within SHERPA are tuned internally, so there are no attributes to define for this method, except the number of evaluations the users wish to perform. SHERPA requires significantly fewer model evaluations than other leading methods do to identify optimized designs, and often finds a solution the first time. This efficiency can save days or even weeks of CPU time during common engineering optimization studies. In the context of calibration, goodness of fit (GoF) is an important measure. It refers to the error between the results of a model and the data that it is trying to replicate (hereafter referred to as “experimental data”). Knudsen and Fotheringham [104] experimented with a number of goodness- of-fit statistics and found the standardized root mean square error (RMSE) to be the best performing. The RMSE can be defined as 2 √(∑(y′ −yi) /n) i RMSE = ̅ (4.1) y where yi′ is the predicted value at matrix point i, yi is the actual value at, ȳ is the mean value of the 85 predicted values (y′), and n is the total number of values. The lower limit of the statistic is 0 which indicates no difference between the predicted values (yi′ ) and the observed values (yi ). The upper limit is usually 1 [104] . RMSE indicates the absolute fit of the model to the data–how close the experimental data points are to the model’s predicted values. Lower values of RMSE indicate better fit. RMSE is a good measure of how accurately the model predicts the response, and it is the most important criterion for fit if the main purpose of the model is prediction. In the current study, an optimization study is conducted for six operating conditions of Cummins Acadia engine with an objective of matching to the experimental second land pressure, third land pressure and blowby using SHERPA method in HEEDS. The RMSE values are observed to find the optimal design at each operating condition. 4.2 Methodology 4.2.1 Problem Setup in CASE The goal of this study is to minimize the RMSE between experimental and CASE predicted the land pressures (second and third) in the piston ring pack of a Cummins Acadia engine. All the data about the geometry, material etc. are known and are used in CASE input file to define the complete cylinder kit. The cylinder-kit contains two compression rings and one oil ring. Objective: • Minimize RMSE of second land and third pressure. Constraints: • Pressure at the end of the entire four stroke cycle should converge by the 9th cycle. Variables: There are ten design variables that are related to the flow paths. The orifice discharge coefficient 86 value is inputted in CASE using the DCOEFF data line. This value is used to define the gas flow equation for gas flow through the ring end gaps. The discharge coefficients for each ring side clearance are defined using the GASFLO input data line. This data line is used to define the gas flow equation for gas flow through the ring side clearances. There are three GASFLO variables for three rings. The equation used for the purpose of both DCOEFF and GASFLO this study was: 𝑚̇ ( 𝐴 ) = 𝐶𝑑 𝜌𝑐 𝜂 (4.2) Where 𝜂 is a factor which is a function of pressure and temperature Cd is the discharge coefficient (DCOEFF and GASFLO) used as a variable in the study. The third type of coefficient is “LEAKAGE” coefficients. Whenever a ring is seated against the top or bottom of its groove, a perfect seal is not always obtained. The LEAKAGE data line is used in CASE to define the sealing coefficients between rings and grooves by defining the leakage area between the ring and its groove whenever it is seated at the top or bottom of its groove due to the imperfect seal and non-regularity of contact surfaces. The leakage area ((𝑃𝑊)𝑙𝑒𝑎𝑘 ) is calculated as following: (𝑃𝑊)𝑙𝑒𝑎𝑘 = 𝐾(𝑖, 𝑗) × 𝑃𝑊 (4.3) where, 𝑃𝑊 is the area between the ring and upper flank of the groove for a bottom seated ring, leakage ratio, 𝐾(𝑖, 𝑗) is the leakage ratio between upper or lower side of the ring and upper or lower side of the groove. In Figure 4.1, the leakage area calculation of a hypothetical ring pack geometry is illustrated. In Figure 4.1(a) all the rings are bottom seated with no leakage. Figure 4.1(b) shows that the rings 87 are lifted which creates leakage area between the ring top and groove top flank and ring bottom and groove bottom flank. Figure 4.1 Ring pack geometry with (a) no leakage in between second ring and the bottom flank of the groove (b) leakage between second ring and the bottom flank of the groove [29] The leakage ratio for the upper and lower side of each ring is defined using the LEAKAGE data line in CASE. The data line allows the user to specify a gas leakage area through the groove when a ring is seated against its groove due to the imperfect seal and non-regularity of contact surfaces. The leakage area is calculated using the leakage ratio as following: (PW)leak = K (i, j) * PW (4.4) where, K(1,1) Leakage ratio for upper side of top ring/groove, K(1,2) Leakage ratio for lower side of top ring/groove K(2,1) Leakage ratio for upper side of second ring/groove K(2,2) Leakage ratio for lower side of second ring/groove K(3,1) Leakage ratio for upper side of third ring/groove K(3,2) Leakage ratio for lower side of third ring/groove K(4,1) Leakage ratio for upper side of oil ring/groove K(4,2) Leakage ratio for lower side of oil ring/groove 88 The values for K(3,1) and K(3,2) are kept as 1 since this cylinder-kit has two compression rings. The rest of the coefficients are tagged as variables in the input DAT file. All the variables used are dimensionless. Since CASE has the option for 3 compression rings and an additional oil ring, the value for each variable for the 3rd compression ring is kept constant as 1 because the piston pack used for this study comprised of only the top, 2nd and the oil ring. Table 4.1 shows the range of variables. Table 4.1 Range for variables Variable Minimum Maximum Orifice discharge coefficient (DCOEF) 0.2 0.7 Discharge coefficient for top compression ring side clearance (GASFLO1) 0.01 0.7 Discharge coefficient for second compression ring side clearance (GASFLO2) 0.01 0.7 Discharge coefficient for oil control ring side clearance (GASFLO4) 0.01 0.7 Leakage ratio for upper side of top ring/groove (LEAKAGE1) 0.0001 0.05 Leakage ratio for lower side of top ring/groove (LEAKAGE2) 0.0001 0.05 Leakage ratio for upper side of second ring/groove (LEAKAGE3) 2e-6 0.01 Leakage ratio for lower side of second ring/groove (LEAKAGE4) 2e-6 0.01 Leakage ratio for upper side of oil ring/groove (LEAKAGE7) 0.0001 0.01 Leakage ratio for lower side of oil ring/groove (LEAKAGE8) 0.0001 0.01 4.2.2 Optimization Study Setup in HEEDS To run the optimization study in HEEDS, the input “DAT” file with all the input data lines used to run a successful CASE simulation is used. The output file used for this study is the “.P02” output file which contains all of the inter ring pressure data for all three rings, and the sump pressure. This output file can be found in the “RPLOT” folder after a simulation is run on CASE. After the 89 input and output files are recognized, an execution command is specified in the execution tab. The pressure convergence output is obtained from an output file called the “OPTIMIZE.RI” file. This file also contains other important data like the average mass flow out of the rings, the total wear per cycle, friction loss, blowby etc. For the purpose of this study, only the pressure convergence is tagged in this output file. The “OPTIMIZE.RI” file can be found in the “RING” folder after a simulation has been run. Since both the output files are inside two different folders, a child folder of the analysis folder is specified in each of the output file properties. After everything related to CASE was set up in HEEDS process tab, the variables and responses are defined according to the problem statement. The resolution for each variable is kept 100001 for accurate results. The number of iterations is set to 1000 for each study. The response for pressure convergence is tagged in the “OPTIMIZE.RI” file and the responses for 2nd land pressure and 3rd land pressure in the “.P02” file. There are a total of 721 output values for both land pressures corresponding to each crank angle of 0 to 720 degrees. To tag the land pressure responses, the script method is used in tagging the entire columns for both 2nd and 3rd land pressures. The next section discusses the land pressures and blowby result comparison between the experimental and CASE predicted results. 90 4.3 Results and Discussion 4.3.1 Optimization at all operating conditions using ten variables. Individual optimization study is conducted for each of the six operating conditions with ten variables. The best design values for all ten variables are listed for each operating condition in Table 4.2 and Table 4.3 shows the resultant blowby. Table 4.2 Best Design values of all ten variables for different operating conditions Operating Best Design Value Condition DCOEF GASFLO LEAKAGE 1 2 4 1 2 3 4 7 8 1800 RPM 0.34 0.18 0.05 0.7 0.05 0.008 2e-6 0.007 0.002 0.00099 100% 1800 RPM 0.36 0.2 0.135 0.55 0.03 0.016 0.0065 0.0098 0.0005 0.0005 51% 1200 RPM 0.27 0.09 0.04 0.1 0.05 0.006 0.0099 0.0094 0.01 0.0001 51% 1200 RPM 0.695 0.02 0.13 0.25 0.0006 0.0002 0.00994 0.006 0.0098 0.002 0% 650 RPM 0.46 0.41 0.07 0.49 0.015 0.002 0.0003 0.006 0.009 0.005 100% 650 RPM 0.41 0.18 0.035 0.47 0.0001 0.005 0.009 8.1 e-6 0.004 0.005 51% At 1800 RPM 100% load, the second land pressure [Figure 4.2 (a)] matches the experimental pressure well. The third land pressure [Figure 4.2 (b)] follows the same trend, but the CASE predicted curve presents a jagged pattern. The blowby is predicted to be 197.05 L/min while the experimental blowby for 1800 RPM 100% is 145 L/min. At 1800 RPM 51% load the second land pressure follows the experimental land pressure nicely, except for a drop and a second peak at 460 CA, see Figure 4.3 (a). Figure 4.3 (b) shows the third land pressure, which presents a jagged pattern. The blowby at this condition is predicted to be 231.15 L/min, the experimental result is not available. 91 Table 4.3 Blowby and gas flow coefficients used across different operating conditions Operating GasFlo_Top GasFlo_2nd GasFlo_Oil Blowby (L/min) Condition 1800_100% 0.18102409 0.0543622 0.7 197.05 1800_51% 0.2 0.135 0.55 231.15 650_100% 0.406591 0.0667663 0.488839 82.74 650_51% 0.412915 0.183356 0.473287 63.68 1200_51% 0.0901504 0.0377656 0.103474 63.75 1200_0% 0.181024 0.0543622 0.7 90.39 Figure 4.2 Comparison between experimental and CASE predicted results at 1800 rpm 100% load (a) 2 nd land pressure (b) 3rd land pressure. 92 Figure 4.3 Comparison between experimental and CASE predicted results at 1800 rpm 51% load (a) 2nd land pressure (b) 3rd land pressure. At 1200 RPM, only the land pressure experimental data at 51% and 0% load conditions are available. No blowby results are available. At 51% load, both the second land pressure [Figure 4.4(a)] and the third land pressure [Figure 4.4 (b)] match the experimental pressure well. At 0% load the second land pressure peak value goes to 520kPa, whereas the experimental peak pressure is 380 kPa, see Figure 4.5(a). Figure 4.5 (b) shows the third land pressure, which presents a reasonable match to the experimental land pressure. The CASE predicted blowby at these two operating conditions are 82.74 L/min and 63.68 L/min, respectively. 93 Figure 4.4 Comparison between experimental and CASE predicted results at 1200 rpm 51% load (a) 2nd land pressure (b) 3rd land pressure. Figure 4.5 Comparison between experimental and CASE predicted results at 1200 rpm 0% load (a) 2 nd land pressure (b) 3rd land pressure. At 650 RPM 100% load, the second land pressure [Figure 4.6 (a)] matches perfectly with the experimental pressure. The third land pressure [Figure 4.6 (b)] does not match. Note that, there the third land experimental pressure presents a constant pressure from 360 CA to 510 CA, indicating experimental discrepancies. The blowby is predicted to be 63.75L/min. At 650 RPM 51% load the second land pressure demonstrates a good match with the experimental results as seen in Figure 4.7 (a). Figure 4.7 (b) again shows the discrepancy in the experimental third land pressure. The 94 blowby at this condition is predicted to be 90.39 L/min, the experimental result is not available. Figure 4.6 Comparison between experimental and CASE predicted results at 650 rpm 100% load (a) 2 nd land pressure (b) 3rd land pressure. Figure 4.7 Comparison between experimental and CASE predicted results at 650 rpm 100% load (a) 2 nd land pressure (b) 3rd land pressure. 4.3.2 Optimization at five operating conditions using the gas flow coefficient used for 1800 rpm 100% The 1800 RPM 100% load is the only operating condition for which the experimental blowby is available. After calibrating that condition the gas flow coefficients for the three rings 0.18102409, 0.0543622 and 0.7 are applied to the rest of the operating conditions. The best design values for 95 seven variables and the resultant blowby are listed for each operating condition in Table 4.4. Table 4.4 Best design values of seven variables while using the same gas flow coefficients as in 1800 RPM 100% in the other operating conditions Operating Best Design Value Blowby (L/min) Condition DCOEF LEAKAGE 1 2 3 4 7 8 1800 0.34 0.05 0.008 2e-6 0.007 0.002 0.00098 197.05 RPM_100% 1800 0.46 0.001 0.0118 0.0034 0.002 0.008 0.0037 184.4 RPM_51% 1200 0.34 0.015 0.0037 0.006 0.00015 0.00598 0.0001 54.96 RPM_51% 1200 0.44 0.0076 0.0069 0.0058 0.01 0.0028 0.005589 67.56 RPM_0% 650 0.54 0.001 0.0028 0.00068 0.002 0.0002 0.0047 89.34 RPM_100% 650 0.42 0.0001 0.0066 0.0015 0.00996 0.0001 0.00698 104.96 RPM_51% Figure 4.8 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 1800 rpm 100%) at 1800 rpm 51% load (a) 2nd land pressure (b) 3rd land pressure At 1800 RPM 51% load, the second land pressure [Figure 4.8 (a)] matches the experimental pressure well. The third land pressure [Figure 4.8 (b)] follows the same trend, but the CASE 96 predicted curve presents a jagged pattern. The blowby is predicted to be 184.4 L/min. Figure 4.9 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 1800 rpm 100%) at 1200 rpm 51% load (a) 2nd land pressure (b) 3rd land pressure At 1200 RPM 51% load, both the second land pressure [Figure 4.9(a)] and the third land pressure [Figure 4.9(b)] match the experimental pressure well. At 0% load the second land pressure drops and peaks at 450 kPa, see Figure 4.10(a). Figure 4.10(b) shows the third land pressure is lower than the experimental land pressure. The CASE predicted blowby at these two operating conditions are 54.96 L/min and 67.56 L/min, respectively. 97 Figure 4.10 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 1800 rpm 100%) at 1200 rpm 0% load (a) 2nd land pressure (b) 3rd land pressure Figure 4.11 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 1800 rpm 100%) at 650 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure At 650 RPM 100% load, the second land pressure [Figure 4.11(a)] matches perfectly with the experimental pressure. The third land pressure [Figure 4.11(b)] does not match and shows the experimental discrepancies. The blowby is predicted to be 89.34 L/min. At 650 RPM 51% load the second land pressure demonstrates a good match with the experimental results as seen in Figure 4.12(a). Figure 4.12(b) again shows the discrepancy in the experimental third land pressure. The blowby at this condition is predicted to be 104.96 L/min, the experimental result is not available. 98 Figure 4.12 Comparison between experimental and CASE predicted results (using gas flow coefficient same as 650 rpm 51%) at 650 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure 4.4 Summary and Conclusions In this chapter, an extensive optimization study is conducted on a range of operating conditions of the Cummins Acadia engine using CASE and HEEDS. The key findings are- • 1800 RPM 100% is the only data set that had experimental blowby results available. • A study is conducted for the 1800 RPM 100% targeting the 2nd land pressure and 3rd land pressure. A very wide range was set for all 10 variables. This study resulted in a reasonable agreement between the CASE predicted and experimental values of blowby, 2nd land pressure and 3rd land pressure. • Using the gas flow coefficients from the 1800 RPM 100% for five other operating conditions resulted in a good match with the experimental second land pressure. Third land pressure prediction followed the trend in all conditions except 650 RPM (both 51% and 100%). The experimental third land pressure at 650 RPM was not reliable. Blowby 99 validation for these five conditions was not possible since experimental blowby was not available. • Obtaining an exact match to the blowby value is possible, but that would compromise the land pressures. • Targeting only the land pressures results in a better overall agreement between the experimental and model prediction. Adding blowby to the list of objectives might give a comparable blowby, but the land pressures would be way off. • No best design is found with a second ring gas flow coefficient of over 0.2 across different operating conditions. 100 Chapter 5 EXPERIMENTAL VALIDATION OF THREE-DIMENSIONAL MULTIPHASE PHYSICS-BASED MODELING METHODOLOGY FOR ENGINE CYLINDER-KIT ASSEMBLY 5.1 Abstract The piston ring pack dynamics and the gas flow dynamics are critical to engine cylinder-kit tribology and design considerations. The current study is the continuation of a three-dimensional multi physics cylinder kit modeling methodology development effort to better understand the performance and physics inside the cylinder-kit. We developed a three-dimensional (3D), multi- phase, multi-physics methodology to investigate the liquid oil- combustion gas transport and oil evaporation mechanisms inside the whole domain of the cylinder kit assembly during the four- stroke cycle. In the current study, the modeling methodology is applied to a Cummins 6-cylinder, 137.02 mm bore, Acadia engine operating at 1800 RPM, full load. First, a CASE (Cylinder-kit Analysis System for Engines) 1D model is developed and a detailed optimization study is conducted to calibrate the model. Next, the ring-bore and ring groove conformability along with the twist angle variation across the circumference are investigated by modeling the positive twisted second ring via a ring FEA contact model. The ring twist induces change in ring location which subsequently changes the cylinder kit geometry dynamically across the cycle. The dynamically varying geometries are generated using a program named LINCC (Linking CASE to CFD) developed in house. Finally, a three-dimensional multiphase flow model is developed for the dynamic geometries across the cycle using CONVERGE. The blowby, second land pressure and third land pressure are compared to the experimental results. The reverse blowby and oil consumption are computed. The liquid oil and oil vapor mass fraction distribution pattern across the cycle are also analyzed. The proposed methodology will elucidate new physics, promote new cylinder-kit designs and has the potential to influence manufacturing of these systems on a grand 101 scale. 5.2 Introduction The present study utilizes the three-dimensional multiphase physics-based modeling methodology presented by Chowdhury et al. in [83] and Part II [105].In Part I [83], the details of the model development methodology have been outlined to explore liquid oil-gas (combustion gases and oil vapor) transport inside the cylinder-kit assembly. Part II [105] presented the effect of a negatively twisted ring on the overall liquid oil-gas (combustion gases and oil vapor) transport along the piston ring pack in a three-dimensional model for the first time. In the current study the methodology is applied on a Cummins 6-cylinder Acadia Engine operating at 1800 RPM, full load. The cylinder kit of this engine has a keystone top ring, a positively twisted taper-shaped ring and a two-piece oil ring. The overall liquid oil-gas (combustion gases and oil vapor) transport along the piston ring pack is explored using the three-dimensional modeling methodology [83,105]. The land pressures and blowby results are compared to the experimental ones. The reverse blowby and oil consumption are also computed. The general flow pattern, liquid oil and oil vapor distribution across the piston ring pack are also analyzed. The entire workflow of the model development consists of four distinct stages outlined as Block 1(One dimensional model in CASE), Block 2 (Three-dimensional Finite Element Analysis (FEA) contact model), Block 3 (Dynamically varying three-dimensional geometry generation), and Block 4 (Three-dimensional Computational Fluid Dynamics (CFD) model in CONVERGE), in Figure 5.1. The inputs, outputs and development details of each block is described in the following subsections. 102 5.3 Workflow Figure 5.1 From tools to the solution (workflow) 5.3.1 Block 1(One dimensional model in CASE) Figure 5.2 Schematic diagram showing three rings with the corresponding ring face profiles In the first stage, following the procedure as described in Part I and Part II [83], a quasi-one- dimensional model is developed in CASE [29], applying a positive static twist of +1.0 to the taper- 103 shaped second compression ring. Primary engine dimensions and operating conditions are presented in Table 5.1. Table 5.1 Engine operating conditions Parameter Value Bore 137.02mm Stroke 169mm Connecting rod length 261.5 mm Engine speed, load 1800 rpm, full load Compression ratio 17.2 Engine oil SAE 10W-30 Cylinder-kit geometrical dimension, material properties, flow characteristics, thermodynamic attributes are some key input parameters for the model development in CASE, as depicted in Block 1 of Figure 5.1. The piston ring pack consists of a low alloy steel keystone top ring, heat treated cast ductile iron taper shaped second ring, cast ductile iron two-piece oil ring and a simplified piston skirt. The schematic of the ring pack with ring running face is shown in Figure 5.2. The in-cylinder pressure and temperature inputs are presented in Figure 5.3 and Figure 5.4 respectively. 104 Figure 5.3 In-cylinder pressure input to CASE model Figure 5.4 In-cylinder temperature input to CASE model The model also considers the liner temperature and distortion, see the liner temperature and 105 distortion in Figure 5.5 and Figure 5.6 respectively. In an additional study, the model is calibrated using HEEDS by varying 10 coefficients such as- 1. Discharge coefficient (1) of three rings 2. Gas flow coefficient (3) for three ring groove side clearance 3. Leakage coefficient (6) for both upper and lower sides of each ring (Modify this with CASE manual) Figure 5.5 Liner temperature 106 Figure 5.6 Cylinder liner deformation The optimization study is conducted in HEEDS. The blowby for one cylinder is found to be 32.84 L/min and accordingly, the blowby for six cylinders is estimated to be 197.05 L/min compared to the experimental blowby value of 145 L/min. The second land and third land pressure curves show good agreement with the respective experimental pressures, see Figure 5.7. The calibration is carried out for five other operating conditions namely- 1800 rpm 50% load, 1200 rpm no load, 1200 rpm 50% load, 650 rpm 100% load and 650 rpm 100% load. The crank angle resolved piston- ring motion from the CASE analysis is exported to be the boundary conditions for the second stage, i.e., 3D ring FEA contact model. 107 Figure 5.7. Comparison between experimental and CASE predicted results at 1800 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure 5.3.2 Block 2 (Three-dimensional Finite Element Analysis (FEA) contact model) In the second stage, a three dimensional FEA contact analysis [39,79] is conducted for the positively twisted second ring. The details of the model inputs are described below. 5.3.2.1 Ring Section Properties As shown in Block 2 of Figure 5.1, the ring cross-section geometry, ring material and thermal properties, inter-ring pressure, ring velocity, and ring acceleration are fed into a three-dimensional ring FEA contact model developed by Cheng et al. [39,79]. The key ring parameters are presented in Table 5.2. 108 Figure 5.8 Temperature boundary conditions of the second ring Under the temperature boundary conditions shown in Figure 5.8, the temperature distribution for the ring cross-section is shown in Figure 5.9. Figure 5.9 Ring cross section temperature distribution Under the current boundary condition, it can be found that the temperature is higher at the ring top outer corner and lower at inner bottom corner. The average temperature is found to be 149.03 °C. As the thermal boundary condition is axisymmetric, the temperature distribution is identical for 109 each cross-section. Table 5.2 Main parameters for the ring Ring material 41114 KV4 -heat treated cast ductile iron according to ISO 6621-3 subclass MC 56 Cylinder bore diameter 137.02 mm Modulus of elasticity 150 GPa Poisson’s ratio 0.3 Density 7150 Kg/m3 Coefficient of thermal expansion 1 E-05 1/K Thermal Conductivity 43 W/mK Ring/gas convective coefficient 25 W/m2 K Ring/oil film convective coefficient 100 W/m2 K Figure 5.10 Second ring motion across the cycle input to 3D FEA contact model (a) ring velocity (b) ring acceleration The ring velocity and acceleration inputs to the contact model are obtained from the CASE model described in block 1, see Figure 5.10(a) and Figure 5.10(b), respectively. 110 5.3.2.2 Experimental Inter-ring Pressure For the current study, the inter-ring pressure is obtained from experiments conducted at Cummins Inc. Note that for the second ring, inter-ring pressure means the pressure at the adjacent boundaries, namely- 1. Second land- This is the surface between the top and second groove 2. Second groove- The taper faced second compression ring is inserted in this groove 3. Third land- This is the surface between the second and third groove Figure 5.11 Inter-ring pressure for second ring input to 3D FEA contact model generated by CASE The pressure across the cycle at these three boundaries with a schematic showing the boundaries is presented in Figure 5.11. The in-cylinder pressure is superimposed for reference. 111 5.3.2.3 Free Shape Ring Measurement Procedure Ring outer diameter (OD) profile at its free state known as ring free shape profile is an input to the model. This is important to accurately predict how the ring will conform to the bore. A piston ring tension test rig [106] is modified to take free ring shape measurements. The initial test rig utilizes 5 micrometer and load cell combination probes to apply various loads to piston rings, see Figure 5.12(a). To measure the free ring shape of the Acadia second ring, additional channels are added to move the probes further out from the center reference. With these new channels added, the test rig can accommodate free ring shapes ranging from 3.35” (85.09 mm) to 6.60” (167.64 mm). These channels are milled with a rotary table to ensure that they remain concentric to the table’s center position. In addition, a new reference puck is machined to a diameter of 137.02 mm corresponding to the Acadia engine’s bore diameter. To measure the free ring points, two probes are held stationary while the third probe is rotated in increments of 5 degrees. This is done on one half of the ring assuming that the ring is symmetrical between halves. Before placing the ring on the rig, the new reference puck is used to zero the probes to the bore diameter. To position the ring, an angled stop (in Figure 5.12(c)) is designed based on measurements of the free shape gap. This angled stop constrains the ring from both rotating and translating. The probe 180 degrees away from the ring gap is left stationary at the bore diameter. The two additional probes are then initially positioned 90 degrees away from the ring gap (180 degrees apart from each other). Then the probes are moved simultaneously in .001 mm increments until their load cells measure a force. At this location, the stationary probe [Figure 5.12 (c)] is locked in as this is the free shape limit for that location. With the ring fully constrained, the moving probe [Figure 5.12(c)] is then moved in 5- degree increments starting furthest away from the ring gap. The micrometer is tightened inward until a force is measured on either of the other two probes, when this occurs, it is immediately 112 stopped, and the location is recorded. This location is relative to the bore diameter since that is what the probe was previously calibrated to and thus is added to that reference diameter of 137.02mm. The calibration of the traveling probe is also checked by first removing the ring and angled stop and then replacing them with the reference puck every 15 degrees. The traveling probe is then re-zeroed if needed before the ring and angled stop are replaced. This swap may introduce some degree of error to the measurement process, but it is required to ensure the micrometer zero is not drifting as it is constantly rotated in and out. Figure 5.12 Test rig for ring free shape measurement. (a) Initial Test Rig for Ring Tension Measurement (b) Test Rig Modifications for Free Ring Shape Measurement (c) Ring Free Shape Measurement Following measurements, the data is input into the Minitab statistical program. Here a fitted line plot is constructed to act as both an averaging mechanism and to complete measurement gaps where the probe could not measure at. The best match to the data was a third-degree polynomial with an R-squared value of 99.9% showing that the regression model accounted for 99.9% of the variance in the data. This is especially useful at the low degrees where the back probe interferes with the traveling probe. The final ring free shape profile is presented in Figure 5.13. 113 Figure 5.13 Second ring free shape profile To validate the reliability of the test rig, a three-probe configuration is used to measure the ring tension and compare the result with the suppliers’ specified tensions for the Acadia ring. The probes are in 90 degree increments to the ring gap as shown in Figure 5.12(a). The back probe (270 degrees in figure 5.12(a)) is left at the bore diameter for the duration of the test. It functions as a stop to ensure the ring is positioned correctly between the remaining two probes. A total of eight trials were conducted with the measurements averaged together. The comparison to the drawing values is shown in table 1. Table 5.3 Acadia Ring Tension Measured tension Supplier provided tension Supplier provided tension Ring (Tangential, N) min (Tangential, N) max (Tangential, N) Top ring 26.73 19.7 29 2nd ring 26.2 17.9 29.8 It is observed that the measured tensions are within the tolerances specified by the ring supplier. This comparison confirms the validity of the developed test rig and the free ring shape measurement conducted using this. 114 The ring dynamics over an engine cycle is analyzed using the finite element method with eight- node hexahedral elements by employing the penalty method [86]. The 3D analytical ring model [39,79] addresses the pressure/force distribution and ring conformability between the ring-cylinder bore interface and ring-groove interface. The contact pairs between ring-liner and ring-groove interface are computed via the penalty method [86] and the force release technique. Note that, bore deformation is not considered in the model. For each crank angle, the twist angle variation over circumferential locations, deformed ring shape profile, piston ring OD profile, and ring location is calculated. The twist angle vs. circumferential location plot at each crank angle is carefully inspected to find out how the ring twist angles are affecting the ring geometry dynamically. Figure 5.14 shows the general trend of the twist angle variation which remains same at each stroke, for one half of the ring. The other half is assumed to be symmetric. The ring twists positively from the ring back to the ring end. Five circumferential locations are picked to characterize the twist angle variation, namely 1, 45, 88, 132, and 168 degrees. These are denoted by A, B, C, D, and E, respectively, in Figure 5.14. Note that 0 degree refers to the ring back, and 168 degree refers to the ring end gap. It is observed that, from the ring back, the twist angle first increases and then drops at 45 degree. It follows the same trend from 45 degree (B) to 88 degree (C). Finally, a sharp drop is observed toward the ring end at 168 degrees (E). 115 Figure 5.14 General trend of twist angle vs. circumferential location After close observation of the variation of the twist angle over the circumference, nineteen (19) crank angles are selected. The twist angle variation is significant at those 19 crank angles. This means the ring location and ring outer diameter (OD) profile undergo a transformation as well throughout the cycle. The selected 19 crank angles at which the ring location and OD change significantly are listed in Table 5.4. 5.3.3 Block 3(Dynamically varying three-dimensional geometry generation) Figure 5.15 Section cut view of the geometry generated using LINCC 116 The ring location and OD profile at the 19 crank angles listed in Table 5.4 are then fed to the next stage LINCC [40], as shown in Block 3 of Figure 5.1 . As mentioned earlier, the geometry changes 19 times throughout the cycle. Those twenty-four dynamically changing geometries are generated using the ring location and OD profiles from the contact model and other geometrical attributes of the cylinder kit. Bore or piston deformation is not considered here. Next, the parts are assembled to form the geometries of the current study in Stereolithographic (STL) format. Figure 5.15 shows a section cut view of one such geometry in Stereolithographic (STL) format. The yellow arrows indicate the generic flow path. The flow domain in the created geometries would be solved in the next stage. Therefore, these dynamically changing geometries from LINCC are used to develop a three-dimensional multiphase model in CONVERGE [46], as shown in Block 4 of Figure 5.1. Table 5.4 Selected crank angles at which the ring location and OD change CA CA CA CA 0 364 489 620 90 374 526 271 389 545 319 418 561 344 455 583 355 472 604 5.3.4 Block 4 (Three-dimensional Computational Fluid Dynamics (CFD) model in CONVERGE The simulation starts at 0 CA with the first geometry and does not change until 181 CA. At 180.99 CA, the results are mapped as the initial condition for the second geometry, and the process goes on throughout the cycle. The quality of the generated mesh is an important parameter in solving the flow physics of internal 117 combustion (IC) Engines. The grid sensitivity study that we conducted in [83] showed that the grids need to be sufficiently fine to resolve the oil-air interface and to capture the tiny ring-groove clearances on the micron scale. Based on that study, the grid size is set at 150 µm in x, y and z- direction, resulting in 3.9 million cells in the computational domain, see Figure 5.16.The simulation was conducted at MSU HPC using 500 CPUs for 1564 hours, resulting in 782,000 CPU hours. It took almost for months to complete the simulation. The mesh is generated automatically in CONVERGE using a cut-cell cartesian meshing technique resulting in a body-fitted volume mesh. The interior cells are hexahedrals. The cells intersected by the geometry surface are cut by the surface to form arbitrary-sided polyhedra [46]. Figure 5.16 Generated mesh of the domain in CONVERGE [46] To analyze the liquid oil-gas multiphase flow, the volume of fluid model (VOF) [48] is used. It is assumed that the gas phase is a combination of N2 and O2 while SAE 10W30 oil is the liquid phase. The flow is assumed to be laminar. The homogeneous relaxation model [57] is used to model the evaporation, assuming that n-heptane (C7 H16 ) is the most volatile constituent of the liquid oil. For liquid oil- combustion gas mass fraction initialization, the domain is divided into 118 four regions such as top ring region, 2nd ring region, oil ring region and skirt region as shown in Figure 5.18. Each ring region consists of the respective land, groove, ring and adjacent part of the liner. The skirt region consists of the piston skirt and the adjacent part of the liner. At the beginning of the cycle, liquid oil mass fraction Yoil is initialized to be 0.05, 0.1, 0.15, and 1 at the top ring region, second ring region, oil ring region, and skirt region, respectively (details in [83]). The temperatures at different regions of the piston and liner are illustrated in Figure 5.17. The liner temperature at the four liner regions are average from Figure 5.5. The in-cylinder pressure and temperature boundary conditions across the domain are shown in Figure 5.3 and Figure 5.4, respectively; see [83] for the detailed description of the methodology including boundary conditions and governing equations. Figure 5.17 Schematic View of the Boundary and Initial Conditions 119 Figure 5.18 Section view of the cylinder-kit with initialized oil mass fraction values 5.4 Results & Discussion 5.4.1 Circumferential flow To compare the general circumferential flow behavior, the stream traces on the three lands of the positively twisted cylinder-kit geometry at a selected CA (Crank Angle) 418 are depicted in Figure 5.19(a), Figure 5.19(b), and Figure 5.19(c), respectively. It is observed that the flow velocity across the lands is influenced by the ring end gap in both untwisted and twisted configuration. In the top land, flow is mostly axial [Figure 5.19(a)]. Figure 5.19 (b) shows the flow in the second land with a view of the top ring end gap and second ring end gap. The flow in the second land is greatly affected by the respective location of the first and second ring end gaps. The flow is axial ranging between 25 m/s-50 m/s just below the top ring end gap whereas the small vortices ranging between 120 25m/s- 62m/s are generated above the second ring end gap. There are few traces with a velocity of 87m/s and above. Figure 5.19 Stream traces at a selected crank angle (418 CA) (a) top land (b) second land (c) third land The flow circulating structures exiting the second ring end gap into the third land are smaller than on the first or second lands. At both sides of the ring end gap the flow structures have circumferential components, but far away from the gap the flow is found to be axial. Note that, due to time constraint, this simulation has been run at a mesh size of 150 microns whereas the nominal piston-bore clearance is 45 microns. Clearly this mesh size is inadequate to capture the actual flow phenomena inside the cylinder-kit. 121 5.4.2 Blow-by To calculate the blowby, mass flowrate at the outflow boundary is plotted in Figure 5.20. In this study, the engine runs at 1800 RPM, and thus it takes 0.033 seconds to complete a cycle. Therefore, 5 1 CA corresponds to 54 millisecond. The discharged mass in kg is found by integrating the area under the mass flow rate versus time graph. The discharged mass is calculated to be 0.0031 kg. This cumulative mass is the sum of the mass of oil and gas (combustion gases and oil vapor). Observing the oil mass fraction at the outflow boundary in Figure 5.21, the averaged combustion gas mass fraction is found to be 0.0.347 at the outlet boundary. Therefore, the discharged combustion gas per cycle is found to be 3.27e-03 kg/s. The final step to calculate blowby is to convert the mass flow rate from kg/s to liter/min. It is assumed that air density is 1.1766 kg/m3 considering air pressure and temperature to be 101325 Pa and 300K, respectively. Finally, the blowby with the positively twisted second ring is calculated to be 169.3 liters/min. The experimental blowby 24.33 liters/min. So, it appears that the 3D model is predicting 3 times higher blowby. 122 Figure 5.20 Mass flow rate at the outflow boundary Figure 5.21 Oil mass fraction at the outflow boundary 123 5.4.3 Estimating Reverse Blow-by & Oil Consumption from the Cylinder Kit Estimation of reverse blowby and oil consumption requires the measurement of the amount of combustion gas (from the inter-ring pack) and liquid oil flowing (from the crankcase through the ring pack) out of the inflow boundary toward the combustion chamber during a complete cycle. These are calculated by summing the mass flow rate that leaves the inflow boundary over the whole cycle. Figure 5.22 shows the mass flow rate at the inflow boundary throughout the cycle. Positive values indicate flow out of the boundary or backflow, which is assumed to be the reverse blowby; negative values indicate flow into the boundary or forward flow directed toward the crankcase. Following the same procedure described in the previous section, the positive discharged mass per cycle (backflow or reverse flow) is calculated to be 6.3892e-04 kg. This cumulative mass is the sum of the mass of oil and mass of gas (combustion gases and oil vapor). the averaged combustion gas and liquid oil mass fraction are found to be 0.9983677 and 1.632235e-03, respectively. Figure 5.22 Mass flow rate at the inflow boundary 124 Thus, the discharge combustion gas out of the boundary or backflow rate is 0.0193 kg/s. Similarly, the liquid oil flowing toward the combustion chamber per cycle is calculated to be 34.76 mg/s .Typical engine oil consumption is 2 mg/s per liter of engine displaced volume [90]. For the 2.49199 L engine operated at 1800 rpm, the ballpark estimate of oil consumption is 4.98mg/s. Thus, the model predicts seven (7) times higher than the typical value. 5.4.4 Land Pressure Comparison Figure 5.23 Comparison between experimental and 3D model predicted results at 1800 rpm 100% load (a) 2nd land pressure (b) 3rd land pressure To compare with the land pressure distribution found from the experiment, the pressures across the second and third land are averaged circumferentially over a cut surface. The 2nd land pressure and third land pressure comparison are demonstrated in Figure 5.23(a) and Figure 5.23(b) respectively. It is observed that they follow the same trend but the pressure value prediction by the model is about 5-6 times higher for both land pressures. Note that during experiment sensors are placed at a single location on the land to measure pressure of the land that might not be representative of all the locations. 125 5.4.4.1 Pressure Circumferential variation Figure 5.24 Circumferential pressure variation predicted by 3D model (a) 2nd land pressure at 424 CA (b) 3rd land pressure at 493 CA To further investigate the pressure behavior on the lands, the circumferential pressure variations on the second land and third land are plotted at a selected CA. Figure 5.24(a) shows that the 2nd land pressure varies from 4050kPa to 4200kPa across the circumference. The corresponding experimental pressure on the 2nd land is 1403 kPa. So, the model predicted pressure is about three times higher at this CA. Figure 5.24(a) shows that the pressure varies from 4050kPa to 4200kPa across the Again, the on the 3rd land the pressure variation across the circumference is between 1470 kPa to 1476kPa. The corresponding experimental pressure is 137.7 kPa which is about eleven times smaller at this CA. 5.4.5 Liquid Oil and Oil Vapor Mass Fraction Distribution Along the Piston Length The results of the liquid oil and oil vapor distribution in the untwisted piston ring pack and cylinder liner at selected crank angles across the computational domain are introduced and discussed in 126 detail in [83] for an untwisted geometry. In [105] the comparison of the liquid oil mass fraction between the negatively twisted and untwisted configuration is presented. In this section, the liquid oil and oil vapor distribution at selected CAs are presented for the positively twisted Acadia cylinder-kit configuration. Figure 5.25 Oil mass fraction on the lands at (a) 0 CA (b) 180 CA (c) 258 CA (d) 361 CA (e) 418 CA (f) 720 CA (g)A 2D sketch of the geometry showing the regions of the computational domain At the beginning of the cycle, liquid oil mass fraction 𝑌𝑜𝑖𝑙 is initialized to be 0.05, 0.1, 0.15, and 1 at the top ring region, second ring region, oil ring region, and skirt region, respectively (details in [83]). Figure 5.25(a) shows the initial distribution along the three lands of the piston at 0 CA. The lands are shown in a 2D sketch of the geometry in Figure 5.25(g). At the end of the intake stroke, i.e., 180 CA [Figure 5.25(b)], the oil transport effect due to ring down scraping and pressure gradient is observed. The third land has 𝑌𝑜𝑖𝑙 between 0.1125 and 0.135 which is lower than the initialized value. At the second land there are oil puddles of 𝑌𝑜𝑖𝑙 =0.09. The oil on the top land is mostly transported to regions below leaving oil mass fraction as low as 0. 127 0225.During the compression and early expansion strokes, it is observed that oil is flowing toward the crankcase due to the blowby effect, see Figure 5.25(c). The top two lands have oil mass fraction 𝑌𝑜𝑖𝑙 = 0.0225 .The oil mass fraction on the third land drops as well indicating the net flow toward the crank case. This movement should persist through the early expansion stroke; but due to lack of continuous supply of oil, by 360 CA there is no oil on the lands. The top land is almost dry while the 2nd and 3rd land has 𝑌𝑜𝑖𝑙 = 0.0225 , see Figure 5.25(d). During the late expansion and early exhaust strokes, the effect of reverse blowby is very minimal. But the oil mass fraction rise from almost zero to 0.0225 on the top land indicates that the little amount of oil left on the second land has transported with the reverse blowby through the ring end gap to the top land. Figure 5.25(e) shows this phenomenon at 418 CA when he second land pressure is 4852.863 kPa while the in-cylinder pressure is 3621.486 kPa. At the end of the cycle, all three lands have 𝑌𝑜𝑖𝑙 = 0.0225 as seen in Figure 5.25(f). 128 Figure 5.26 Oil vapor mass fraction on the lands at (a) 180 CA (b) 258 CA (c) 361 CA (d) 418 CA(e)525 CA (f) 720 CA (g)A 2D sketch of the geometry showing the regions of the computational domain Oil evaporation across the piston lands throughout the cycle is illustrated in Figures 5.26(a) through Figure 5.26(f), and the regions are identified in a two-dimensional sketch in Figure 5.26(g). It is observed from Figure 5.26(a) that there is quite some oil vapor on all three lands at the end of the intake stroke. After that, oil evaporation is not prominent during the compression stroke with only small puddles on the second and third land; refer to Figure 5.25(b). At the end of compression stroke there is barely any oil vapor, see Figure 5.25(c), The mass fraction of oil evaporation up to this point is on the order of 10−13 . During expansion and exhaust stroke, oil is transported as vapor into the second land and subsequently to the top land as a result of reverse blow-by, A gradual increase of oil vapor puddles with higher mass fraction on the order of 10−11 is observed in Figure 5.25(d) through Figure 5.25(f). The oil vapor in the top land could be transported to the combustion chamber. 129 5.4.6 Speed Comparison: MSU HPC Vs. Argonne The current simulation has been run on MSU HPC using 500 CPUs for 1564 hours resulting in a total of 782,000 CPU hours. A strong scaling study is conducted in which the problem size n (say, number of grid points) is fixed and the number of cores Pc is increased [107]. The number of cells is kept at 3.9 M with 150 µm grid size. The number of cores is increased gradually from 64 to 512 for both MSU HPC and Argonne THETA. At MSU HPC each user can use a maximum of 520 cores. So, the number of cores is further increased to 8192 at Argonne THETA. Each of the experiment is conducted for one hour. Table 5.5 Simulation completion and network latency comparison for a one-hour clock time: MSU HPC Vs. Argonne THETA Total no MSU: CA MSU: Average Argonne: CA Argonne: Average of CPUs completed computed network completed computed network latency (in latency (in microseconds) microseconds) 64 0.133 3.095324 0.032 10.86 128 0.191 8.933805 0.066 11.78 192 0.324 3.156953 0.107 8.826571 256 0.35 3.144804 0.122 11.157 320 0.378 3.032791 0.138 11.37725 384 0.379 11.974364 0.146 10.85 448 0.372 3.564262 0.17 10.834 512 0.35 12.987704 0.165 11.256 1024 - - 0.222 9.862787 2048 - - 0.257 10.034252 3072 - - 0.23 10.344646 4096 - - 0.202 10.337646 5120 - - 0.194 9.803858 8192 - - 0.135 12.190331 130 Figure 5.27 CAs completed vs. no. of cores during a 1-hour clock time at MSU HPC and Argonne THETA Table 5.5 shows how many crank angles (CA) can be completed in an hour as the number of CPUs are increased at both facilities while everything else of the simulation is kept same. The average network latency at both facilities is also listed in Table 5.5. Figure 5.27 shows at MSU HPC as the number if cores increase more CAs are completed up to 384 cores. After that it drops. On the other hand, the Argonne THETA demonstrates an upward trend as the number of cores is increased. It is also observed that MSU HPC is 3-4 times faster. THETA uses “Intel Xeon Phi x200, codename Knights Landing (KNL), a second-generation MIC architecture product from Intel” and MSU HPC uses “Intel Xeon @ 2.4 GHz+ and AMD EPYC @2.5 GHz+”. A single core of one of these Xeons can be 5-6x faster than a single KNL core. So, we need to use many more KNLs to exceed the performance of a collection of Xeons. 131 Figure 5.28 Average computed network latency vs. no. of cores during a 1-hour clock time at MSU HPC and Argonne THETA The average computed latency at both facilities is also plotted. The network latency indicates how long does it take for the nodes to communicate or transfer data between each other in a network. Since latency varies with trial, the latency numbers are spread in an erratically. It is clearly observed that the minimum and maximum latency of the MSU network are 3 and 12 microseconds, respectively. On the other hand, the minimum and maximum latency of the Argonne network are about 10 and 13 microseconds, respectively. So MSU network would have advantage over Argonne in terms of network latency. Since at Argonne THETA, more CPUs are available, the number of cores is further increased up to 8192. It is observed that THETA scales out after 2000 cores, see Figure 5.29. The network latency also follows the same trend as observed in Figure 5.30. 132 Figure 5.29 CAs completed no. of cores during a 1-hour clock time at Argonne THETA Figure 5.30 Average computed network latency vs. no. of cores during a 1-hour clock time at Argonne THETA 133 5.5 Summary and Conclusions In this study, the three-dimensional multiphysics methodology is applied on the Cummins Acadia engine. A 1D model is developed, using CASE to generate the boundary conditions of the three- dimensional model. The ring free shape is measured using an inhouse test rig and 3D FEA contact model is used to find the crank angles at which the twist angles and ring locations change significantly. Then the dynamically varying geometries at those crank angles are generated using the LINCC program. Finally, a three-dimensional, multi-phase, CFD-model is developed using CONVERGE. • The blowby prediction is about 7 times higher than the experimental blowby. • The land pressures are predicted about 5 times higher than the experimental results. • The one-hour runtime experiment shows that adding number of CPUs is not speeding up the simulation. MSU HPC is observed to perform better for the specific problem. More on the limitations and future recommendations will be discussed in Chapter 7. 134 Chapter 6 TRIBOLOGICAL PERFORMANCE ASSESSMENT OF ABRADABLE POWDER COATED PISTONS CONSIDERING PISTON SKIRT GEOMETRY AND SURFACE TOPOGRAPHY 6.1 Abstract Surface coatings are one of the most widely used routes to enhance the tribological properties of cylinder kits due to effective sealing capability with low friction coefficient and high wear resistance. In the current study, we have conducted the surface texture characterization of the coating on piston skirts and evaluated the impact of a novel Abradable Powder Coating (APC) on cylinder-kit performance in comparison to stock pistons. The surface texture and characteristic properties varying across the piston skirt are obtained and analyzed via a 3D optical profiler and OmniSurf3D software. The engine operating conditions are found through a combination of measurements, testing, and a calibrated GT-Power model. The variable surface properties along with other dimensions, thermodynamic attributes, flow characteristics and material properties are used to build a model in CASE (Cylinder-kit Assembly System for Engines)- PISTON for both an APC coated piston and a stock piston. The surface texture analysis shows that the APC coating resembles a group of mushroom cap structures. The APC coated surface has a higher number of deeper valleys than the stock coated surface; these valleys would potentially be beneficial for lubrication and oil retention. Comparison of different performance parameters from CASE simulation results shows higher minimum oil film thickness and lower hydrodynamic shear force on the APC coated piston. The APC coated piston is predicted to experience smaller tilting as well as secondary motion than the stock piston. All these features are believed to make the APC coating a suitable candidate for piston skirt coating. 135 6.2 Introduction Engine cylinder kit tribology is a very complex phenomenon. Friction, wear, and lubrication processes inside the cylinder kit assembly affect each other through interacting causes and effects, making the issue complex. Interactions between the piston-liner, ring-liner, ring-groove surfaces are highly complex and nonlinear; hence full understanding requires knowledge of friction, wear, and lubrication. Surface coating is one of the most widely used routes to tailor surface morphology and enhance the tribological properties of engine cylinder kit. The coating improves the life of the engine as well as fuel efficiency. Researchers have sorted different surface coatings on piston ring–cylinder liner and piston skirt– cylinder liner interfaces. Experimental studies have focused on investigating the effect of piston ring coating [108–116] and cylinder liner coating [117,118] on friction and wear behavior at their interfaces via simulated laboratory tests. Piston skirt-cylinder liner interface is amongst the major contributors to frictional power losses of Internal Combustion Engines (ICEs). skirt coatings can contribute to minimizing wear and friction by changing the surface topography and controlling clearance and skirt geometry. Federal-Mogul [119] developed a piston skirt coating known as EcoTough-D coating (EcoTough-Diesel). This polymer-based coating was reinforced with short carbon fibers and graphite was embedded in the coating as a solid lubricant. Friction losses were reduced by up to 35 % compared to standard coatings on the market. Demas et al. [120,121] tribologically tested piston skirt segments extracted from a commercial aluminium alloy piston against commercial grey cast iron liner segments using a reciprocating laboratory test rig. The piston skirt segments were coated with diamond-like carbon (DLC) coating, graphite–resin coating nickel–polytetrafluoroethylene (Ni–PTFE), and an amorphous hydrogenated carbon (a-C: H) coating. The friction dependence of these piston skirt 136 and cylinder liner materials was studied as a function of load, sliding speed, and temperature. Wang and colleagues [122] applied nickel–ceramic composites via conventional electroplating over piston skirts and slid against aluminium and cast iron bores to study wear and scuffing characteristics. Gavrilova et al. [123] experimentally compared the friction coefficient of anti- friction solid lubricating coatings dedicated to decreasing friction and wear in the piston skirt– cylinder liner of combustion engines. Solid lubricating coatings based on a polymer binder and high-purity molybdenum disulfide (MoS2) and graphite were chosen for experimental studies. Westerfield et al. [124] discussed friction trends associated with the piston skirt profile as well as coating pattern and roughness of the Mahle GRAFAL coated piston skirt using a floating liner engine experimental study. Cho and colleagues [125] experimentally compared the performance of graphite and diamond-like carbon (DLC) coated piston skirts in terms of friction losses and the amount of wear in a piston assembly. Numerical modeling and simulation can be a time-saving and cost-effective alternative to engine and laboratory testing for the performance assessment of cylinder-kit coatings. Previous works involving modeling and simulation of coated surfaces have studied both ring-cylinder, piston skirt- cylinder interfaces. Zabala et al. [126] focused on coupling the experimental simulation of a ring- on-liner contact at the top dead center (TDC) with computational modeling of wear and friction as a function of chromium-nickel (CrN) physical vapor deposition (PVD) coatings as well as self- lubricious diamond-like carbon and molybdenum disulfide (MoS2) coatings. Wan et al. [116] investigated the interacting mechanism of corrosion and fatigue wear of the chromium-nickel (CrN) coated piston rings by simulating the stress distribution and crack origins in the CrN coating using a finite element method (FEM). Buyukkaya et al. [127] developed a finite element model to conduct thermal analysis on uncoated and ceramic coated diesel pistons, made of aluminum-silicon 137 alloy (AlSi) and steel. Several studies investigated the effect of piston skirt coating on tribological properties via numerical modeling. Hildyard et al. [128] applied three different variants of Graphene Oxide (GO) coatings on aluminium piston skirt of gasoline engine via the electro-phoretic deposition (EPD) method. Their tribological performance was benchmarked against uncoated steel and graphite- coated aluminium skirts. These coatings were experimentally characterized in terms of asperity level friction, topography, and wear resistance. They used those parameters to numerically simulate the performance using a multi-physics tribo-dynamic model developed in AVL Excite Power. Ricardo, Inc. [129] performed simulations to estimate friction reduction and potential fuel economy improvements due to the application of advanced surface treatments at key interfaces of engines such as piston skirt-liner, rings-liner, interfaces within journal bearings, and valvetrain. The Ricardo software PISDYN [130] was used to simulate the friction behavior of the piston skirt as it travels up and down the cylinder liner. Zhang et al. [131] investigated the role of surface and interface morphology of piston skirt graphite coatings on decreasing friction and resisting the scuffing of the cylinder bore. They also carried out series of gasoline tests with different skirt profiles to simulate the friction and wear performance of the skirt with mixed coatings of graphite and MoS2 [132]. In both studies, the dynamic software PISDYN from Ricardo was used to simulate and predict the motion, contact, wear, and kinetic energy in the working process of the piston. Meng et al. [133] modified an existing model of piston skirt lubrication with consideration of breaking in process for the most commonly used triangle machine marks on GRAFAL coated piston skirts. A new set of flow factors in the averaged Reynold’s equation were analytically derived for the trapezoid shape formed after wear of the original triangle shape. A new asperity contact model was developed for the trapezoid shape. 138 Thermal Spray Abradable Coatings have been used for decades to improve the efficiency of engines. Suman et al. [27,28] reported a novel technology of Abradable Powder Coatings™ (APC™) for the application in piston skirts. The materials are based on thermoset powder coating resins, which serve as a binder phase. Resin modifiers and lubricious fillers, such as graphite are added to tailor the properties. The coatings can be applied in thicknesses from 15 µ to over 250 µ as sprayed. The coating is thick and abradable, allowing it to form to its mating surface over time. Generally, a line-to-line fit upon assembly of a device is achieved. During initial operation, the coating breaks in to form a perfect fit between mating parts. Once the optimum fit is achieved, stresses on the coating are reduced and the break-in process stops. The coatings reduce operating clearances on piston, support and maintain hydrodynamic oil film regime. Sliding wear tests and field engine tests showed that this has improved durability compared to liquid piston coatings. Thin, solid film lubricant coatings help after the oil film fails, whereas APC prevent the oil film from failing [28] . APC can take advantage of properties of almost any material by incorporating them into the coating systems as a distribution or as a layer. Key advantages of APC are listed below- 1. APC creates an ideal fit where peak contact stresses are transferred to adjacent areas – lowering the peak stress, and broadening the ‘footprint’ of the load bearing oil film area. A thinner, wider oil film area is more stable and has less turbulent losses than a peaky contact pattern that readily pierces the oil film. 2. On a microscopic level, the plateaued surface of APC coatings is ideal with low peak height and high valley depth providing a polished bearing surface plateaus with surplus oil volume in adjacent pores. 139 3.APC creates a better fit between the piston and bore over time which limits secondary motion of the piston, which is responsible for NVH (noise, vibrations, harmonics), additional wear modes, and parasitic energy losses. 4. The thicker coating of an APC can imbed the foreign particles (from the combustion chamber, bearing wear, casting sand particles, machining burrs, and general dirt that may be captured in the build) into itself. As a result, the particles cannot penetrate the oil film as easily whereas thinner coatings do not have that advantage. In this study, we analyzed the unique surface topography of worn APC coated piston skirts and compared it to the worn stock piston skirts of a Cummins 2.8L, four cylinder, turbo engine. The specific APC formulation used for the “coated” piston skirt in this study is known as “Additive Abradable Graphite Coatings (AAGC)”. The stock piston skirt also has a porous powder coating. Since in this study our focus is to explore the unique surface characteristics of APC, we are denoting the APC piston skirt to be “coated”, and the piston skirt originally installed in the Cummins 2.8L engine to be “stock” in our discussion. First, the topography of the coatings is analyzed using a three-dimensional (3D) optical profiler and a surface texture analysis software OmniSurf3D [136]. A GT-Power model is developed and calibrated using experimental data to obtain engine operating conditions. Finally, a numerical model is developed in the CASE-PISTON module [29] by Mid-Michigan Research for both coated and stock pistons. The variable surface properties varying across the piston skirt surfaces are implemented in Greenwood-Tripp [43] model inside CASE-PISTON. The fluid film lubrication is analyzed by the solution of the Reynolds equation [41]. 6.3 Piston Surface Morphology Analysis Surface topography plays a significant role in understanding the tribological characteristics such 140 as friction, wear, and lubrication between the mating parts [137–139]. In this study, we focused on the surface topography of the “coated” and “stock” piston skirt surfaces of the cylinder kit. Primarily fifteen surfaces were selected on the stock piston for surface measurements. These are the upper and lower side surfaces of each groove and the upper side, lower side, and face of each ring, as shown in Figure 1. These fifteen surfaces (six groove surfaces and nine ring surfaces) have the stock surface characteristics and it is assumed that those surface properties are constant throughout these surfaces. The cylinder surface properties are obtained from available supplier data and assumed to be constant throughout the surface. The surface properties of the surfaces mentioned are used while developing the model in CASE and not discussed in this paper in detail. Two piston skirts (surfaces 16 and 17 in Figure 6.1) of the Cummins 2.8L engine were coated with APC and two were kept at their original coating or “stock”. After running the engines for 1.5 hours one worn stock piston and one worn coated piston were taken out for surface measurements, see Figure 6.2. All surface measurements and subsequent image analysis followed the same protocol. 141 Figure 6.1 Selected surfaces for measurement Figure 6.2 (a) Worn coated piston (b) Worn stock piston Figure 6.2(a) shows the worn coated piston labeled as “E1P2 Post-run Coated” and Figure 6.2(b) shows the worn stock piston labeled as “E1P3-Post run Stock”. On both piston skirts two 142 circumferential locations A (30 degrees clockwise from the center of thrust and B (center of thrust) are selected. Five evenly spaced axial locations are selected between the spacings below each letter, corresponding to each circumferential location. The axial locations are 16.45 mm, 18.79 mm,23.46 mm, and 25.8 mm away from the piston skirt top. The locations are named A1-A5 and B1-B5 for coated and C1-C5 and D1-D5 for stock piston skirts. Table 6.1 shows the complete list of measurement locations on the piston skirts, see those locations in Figure 6.2. It is also assumed that the properties are the same on the anti-trust side. So, surfaces 16 and 17 in Figure 6.1 are considered to have the same surface properties. Table 6.1 List of the piston skirt measurement locations Location Label Circumferential location Axial location (distance from piston skirt top in mm) no. (Degrees from the center of thrust) Location 1 A1(Coated) 30 16.45 C1 (Stock) B1(Coated) 0 D1(Stock) Location 2 A2(Coated) 30 18.79 C2 (Stock) B2(Coated) 0 D2(Stock) Location 3 A3(Coated) 30 21.1256 C3 (Stock) B3(Coated) 0 D3(Stock) Location 4 A4(Coated) 30 23.4624 C4 (Stock) B4(Coated) 0 D4(Stock) Location 5 A5(Coated) 30 25.7992 C5 (Stock) B5(Coated) 0 D5(Stock) The measurements are taken at Michigan Metrology LLC using NPFlex 3D Optical Profiler (Bruker Corporation), see Figure 6.3. The NPFlex is a microscope in which each objective lens 143 contains a specially designed interferometer. The microscope portion of the instrument provides the “image” of the surface whereas the interferometric portion provides the height information comprising the surface. The primary mode of operation is the vertical scanning interferometer (VSI) technique [140]. Figure 6.3 NPFlex 3D Optical Profiler (Bruker Corporation) [140] Table 6.2 Measurement Protocol for the Optical Profilometry Measurement Attribute Nominal Value Magnification 5.5X Measurement Array Size 640 x 480 Lateral Sampling 1.8 um Field of View 1.2 mm x 0.9 mm Height Resolution < 6 nm 3D (Areal) filter Long Wave Pass, S-Filter = 12 μm Terms Removal Best fit Tilt/Curvature Stylus X ls/lc 8 μm / 0.80 mm Stylus Y ls/lc 8 μm / 0.80 mm Table 6.2 shows the measurement protocol used to measure the piston skirts. For typical machine parts, the standard protocol for the field of view is 1 mm x 1 mm [140]. Following the protocol, the field of view is set to 1.2 mm x 0.9 mm. The long wave pass, S-filter for the 3D measurement 1 is set to be 100 of 1.2 mm, i.e., 12 μm. This filtering is used to ensure that there is no noise 144 narrower than 12 μm wavelength. Figure 6.4 Orientation of the sample for measurement of piston skirt surface For the piston skirt surface measurement, the X-axis is toward the top of the piston, and Y-axis is toward the circumference of the piston as shown in Figure 6.4. In the orientation image of Figure 6.4, the dark lines are the grooves. The images taken are analyzed using the OmniSurf3D software by Digital Metrology [136]. The analysis setting parameters are listed in Table 6.3. The surface might be misaligned due to the placement of the surface under the instrument during measurement. In addition, the measured data may contain underlying geometry that is not part of the texture. OmniSurf3D provides tools for establishing a reference geometry that is removed for texture analysis [136]. In the current study, a fourth-order polynomial is selected as reference geometry. To fill the missing data “Bicubic Fill” is used. Once our surface data has been adjusted for geometry, we need to apply filtering to isolate the wavelengths of interest. We have used a Gaussian short-wavelength filter with a 12 microns cut-off. 145 Table 6.3 Image Analysis Protocol for Omnisurf3D Analysis Attribute Analysis Setting Missing Data Filling Bi Cubic Fill Reference geometry Fourth-order polynomial Short Filter Type Gaussian Short Filter Cutoff 12µm Long Filter Type 2nd Order Gaussian Long Filter Cutoff 800 µm X, Y Spacing 1.8 µm, 1.8 µm The three properties considered to be CASE input are: root mean square (Sq) (RMS) roughness (μm), summit mean radius (SSMR) (μm), and summit density (SSDN). These measurements are taken relative to the reference plane. A reference plane is defined by calculating the plane as the mean of all the points. The root mean square (Sq) (RMS) roughness (μm) is the standard deviation of the data point distances to the reference plane (ISO 25178-2). Another term used in tribological modeling for lubrication and contact behavior is the surface summit (local peak). A "summit" is considered as a local peak that is above the reference plane and is higher than all of its 8 nearest neighbors [136]. The summit mean radius (SSMR) (μm) of each summit is calculated as the average of the radii that passes through the considered summit in the two orthogonal directions X and Y. The mean radius of the summits (SSMR) is considered the average of all summit radii [141]. The summit density (SSDN) is calculated by computing the ratio between the number of summits and the measured area. The inverse of this property is known to be the average area represented by a summit (1/SSDN) (mm2 ), which is the input for CASE. To understand the surface texture better, we would also investigate the “valley volume per unit area (SVO) (μm)” that has not been implemented in CASE yet. If we visualize the surface as a series of extreme peaks, kernel or core, and then extreme valleys, SVO is the volume of the extreme valleys. This parameter is normalized by area, hence the unit is μm. This would be the equivalent 146 of the volume of fluid the valley would hold. The following sections contain a detailed analysis of the surfaces of the worn stock and worn coated piston skirts. 6.4 Comparing Surface Properties of Worn Stock and Worn Coated (Post-run) To compare and evaluate the overall surface texture of the two different surfaces some key parameters are selected. They are the summit radius (SSMR), summit density (SSDN), root mean square (RMS) surface roughness, and valley volume per unit area (SVO). Figure 6.5 The surface topology of the worn piston skirt at side location (a) Coated location1 (A1) (b) Stock location 1 (C1) (c) Coated location 4 (A4) (d) Stock location 4 (C4) (e) Worn coated piston, labeled (f) Worn stock piston, labeled. This section compares the surface texture and these parameters of the worn coated and worn stock piston skirts at the side location (A for coated and C for stock) and at the center location (B for coated and D for stock). 147 6.4.1 Worn Coated vs Worn Stock at Side Location At the piston skirt side location, the piston skirt has little wear. Figure 6.5 presents the zoomed-in view of the worn coated and worn stock piston skirts at two representative axial locations-location 1(16.45 mm below the skirt top) and location 4 (23.46 mm below the skirt top). For the full 3D contour images of two representative surfaces, see Appendix A-1 (coated location 1), A-3 (stock location 1). It is observed from Figure 6.5(a) and Figure 6.5(c) that the coated surfaces have a good number of deep pockets or valleys to hold oil. On the other hand, the stock surfaces have barely any such pockets, see Figure 6.5(b) and Figure 6.5(d). The coated surfaces look more like mushroom caps whereas the stock surfaces contain finer feature. Table 6.4 Comparison of surface properties at side locations for stock (C) and coated (A) piston skirts Axial SSDN SSDN SSMR SSMR Sq Sq SVO SVO locations (1/𝐦𝐦𝟐 ) (1/𝐦𝐦𝟐 ) (µm) (µm) (µm) (µm) (µm) (µm) Type Ctd Stk Ctd Stk Ctd Stk Ctd Stk Location 1 1209 1567 31.825 26.126 6.868 2.819 0.878 0.099 Location 2 1147 1801 30.882 33.955 7.158 2.817 0.808 0.13 Location 3 1034 1512 30.701 26.582 7.694 2.999 0.856 0.107 Location 4 912 1523 24.672 26.114 7.757 3.044 0.731 0.117 Location 5 861 1496 20.722 27.075 7.979 2.998 0.768 0.1 Table 6.4 shows the comparison of surface properties at side locations for stock (C) and coated (A) piston skirts. Here, Sq is RMS Surface roughness (μm), SSDN is summit density (1/mm2 ), SSMR is average summit radius (μm), Smr is material ratio (%), SVO is volume of extreme valleys (μm), Ctd is coated and Stk is stock. Table 6.4 shows that the RMS surface roughness of the coated surface at each location is about 148 three times higher than the stock ones. This means the coated surface will hold more oil. The summit density is about 1.5 times lower at every location for the coated skirt than for the stock skirt. The summit radius for both surfaces is comparable since this location is little worn. Further investigation shows that the valley structures play a significant role in the coated surface in terms of oil retention and lubrication. It is observed that the coated surface has deeper valleys (Figure 6.6 (a)) than the stock surface (Figure 6.6(b)). From Table 6.4 it is observed that the valley volume of the coated surface is about nine times higher than the stock surface valley volume. Figure 6.6 Comparing the Valleys at Side Locations(a) Valleys at Worn Coated Location 1 (A1) b) Valleys at Worn Stock Location 1 (C1) 6.4.2 Worn Coated vs Worn Stock at Center Location At the center location, both surfaces wear away and materials are removed. Figure 6.7 presents the zoomed-in view of the worn coated and worn stock piston skirt at two representative axial locations-location 1(16.45 mm below the skirt top) and location 4 (23.46 mm below the skirt top). For the full 3D contour images of two representative surfaces, see Appendix A-2 (coated location 1), A-4 (stock location 1). In a wear process, the surface continually changes. It starts with peaks of a smaller radius. As wear progresses, it leaves a plateau-like surface with a larger radius. It is 149 observed from Figure 6.7 that, both coated and stock surfaces appear flatter than at the side location shown in Figure 6.6. Again, more pockets to hold oil are observed in the coated surfaces [Figure 6.7(a), Figure 6.7(c)] as opposed to finer features observed in the stock surfaces in Figure 6.7(b) and Figure 6.7(d). Table 6.5 shows the comparison of surface properties at side locations for stock (C) and coated (A) piston skirts. Here, Sq is RMS Surface roughness (μm), SSDN is summit density (1/mm2 ), SSMR is average summit radius (μm), SVO is the volume of extreme valleys (μm), Ctd is coated and Stk is stock. Table 6.5 Comparing Center Locations for Stock (D) and Coated (B) Piston Skirts Axial SSDN SSDN SSMR SSMR Sq Sq SVO SVO locations (1/mm2) (1/mm2) (µm) (µm) (µm) (µm) (µm) (µm) Type Ctd Stk Ctd Stk Ctd Stk Ctd Stk Location 2 4250 3229 75.812 45.803 3.233 1.555 0.936 0.417 Location 3 4225 3607 78.473 49.621 3.632 1.363 1.08 0.356 Location 4 4340 2375 75.636 41.886 3.609 2.231 1.14 0.372 Location 5 4178 2309 73.964 46.546 3.808 2.175 1.057 0.285 Table 6.5 shows that the RMS surface roughness, summit density, and surface roughness of the coated surface at each location is higher for the coated surface than the stock ones. It is observed that the coated surface has deeper valleys [Figure 6.8 (a)] than the stock surface [Figure 6.8(b)] at the center location as well. From Table 6.5 it is observed that the valley volume of the coated surface is about three times higher than the stock surface valley volume. The surface properties such as summit density, summit radius, and surface roughness varying along piston skirt axial and circumferential locations are implemented in the CASE (Cylinder-kit 150 Analysis System for Engines) program to investigate how these properties are affecting the performance of stock and coated pistons. Figure 6.7 Surface topology of the worn piston skirt at center location (a) Coated location 1 (B1) (b) Stock location1 (D1) (c) Coated location 4 (B4) (d) Stock location 4 (D4) (e) Worn coated piston, labeled (f) Worn stock piston, labeled. Figure 6.8 Comparing the valleys at center locations (a) Valleys at worn coated location 1 (B1) b) Valleys at worn stock location 1 (D1) 151 6.5 CASE Model Setup Figure 6.9 CASE model setup workflow The entire workflow of the model setup in Cylinder-kit Assembly System for Engines (CASE) consists of two stages labeled as Block 1(Input parameters and files generation) and Block 2 (Running model using CASE-PISTON), in Figure 6.9. 6.5.1 Block 1(Input parameters and files generation) This study is conducted on a Cummins 2.8L turbocharged four-cylinder diesel engine running at 2000 rpm, 12 bars for about 1.5 hours. Certain engine operating conditions are presented in Table 6.6. 152 Table 6.6 Engine operating conditions Parameter Value Nominal bore diameter 94 mm Nominal piston diameter 93.0656 mm Stroke 100 mm Connecting rod length 158 mm Engine speed, load 2000 rpm, 12 bar Number of valves 4 Rated HP 210 HP @ 3600 RPM Rated torque 385 lb-ft Torque @ 1800 Compression ratio 16.9 Engine oil SAE 15W-40 Piston height 75.946 mm Skirt height 46.228 mm The current study employs GT-Power by Gamma Technologies (GT) as the modeling platform to develop a one-dimensional (1D) engine model of the Cummins four-cylinder turbocharged engine which is calibrated using experimental data. The developed model includes intake/exhaust runners and ports, intake/exhaust valves, direct injection (DI) injectors, intercooler, exhaust gas recirculation (EGR) cooler, and controller circuit, compressor, turbocharger, and wastegate controller. The model is developed by using all the piping geometry collected through direct measurement of the actual setup and the provided turbocharger and compressor maps. Appendix A-5 shows the GT map of Cummins 2.8L turbocharged four-cylinder diesel engine. The solution is based on one-dimensional(1D) fluid dynamics, representing the flow and heat transfer in the piping and other flow components of the engine system. The flow model specifically involves the solution of 1D Navier-Stokes equations [142], namely the conservation of mass, momentum, and energy equations. These equations are solved in one dimension which means all quantities are averaged across the flow direction. For the current simulation, the ‘explicit’ time integration method has been employed. Heat transfer in the engine cylinder is modeled using the WoschniGT heat transfer model, which closely emulates the classical Woshni correlation without 153 swirl. The present study employs Gamma Technologies’ semi-predictive direct injection combustion model. The injection is also controlled through a proportional–integral–derivative (PID) controller that attempts to match engine load based on the target input. Turbine and compressor performance are modeled in GT-SUITE using performance maps from Cummins. To check the overall robustness and accuracy of the model, the model is simulated at over 170 cases ranging from idle to maximum RPM and from base to peak engine load. This is generated by taking the published torque curve by Cummins and back calculating the respective brake mean effective pressure (BMEP). This case sweep featured zero failed solutions and showed a great correlation in achieving the peak power at each RPM point compared to the published values. Later the model was calibrated using the experimental data from dyno testing. The in-cylinder combustion gas pressure and temperature profiles obtained are presented in Figure 6.10 and Figure 6.11, respectively. To obtain cylinder and piston temperature inputs, a thermal finite element analysis is conducted using GT’s Spaceclaim and GEM software. The process begins with inputting the model into GT’s Spaceclaim software where the origin and position of the piston are fully defined relative to the engine block and other engine components. This is then input into the GEM software where the FE mesh is constructed for the CAD model and surface interactions are defined. The minimum surface definitions required are ring interactions, combustion gas interaction, potential cylinder contact areas, and one oil interaction zone. These zones were all selected for the piston using the GT manual as a reference to ensure that the correct surfaces were selected for the calculations on each surface. The internal oil gallery was created by implementing unique oil temperature zones throughout the gallery and the ports. These zones were assigned unique heat transfer coefficients (HTC) by summarizing the HTCs found in [143]. A mesh study was also conducted to find the 154 minimum resolution needed for accurate thermal solutions. This was found to be 5mm and was further used as the inputs to CASE simulations. The resultant cylinder temperature and piston temperature for CASE inputs are shown in Figure 6.12 and Figure 6.13, respectively. The cylinder kit assembly in this study is comprised of two compression rings: a keystone top ring, a tapered front face second ring, and an oil control ring. The geometrical dimensions such as detailed measurements of the piston, rings, and liner, including end clearance of piston ring pack, location of ring end gaps are measured. The materials of different components of the cylinder-kit are listed in Table 6.7. The properties of the listed materials are obtained from literature and mechanical engineering handbooks. Table 6.7 Materials of different components Part Material Cylinder liner Grey Cast Iron Piston Hypereutectic Cast AL Alloy (M126/M138) Top ring SAE 9254, Steel 2nd ring F14 - Alloyed tempered flake cast iron Oil ring SAE 9254, Steel Figure 6.10 Combustion gas pressure 155 Figure 6.11 Combustion gas temperature Figure 6.12 Liner temperature 156 Figure 6.13 Piston temperature Figure 6.14 Piston skirt profile of worn stock and worn coated pistons. (a) vertical trace (b) circumferential trace The vertical trace and circumferential traces of both worn stock and worn coated piston skirts are also measured and presented in Figure 6.14(a) and Figure 6.14(b) respectively. Note that, this does not show diameter or waviness. All measurements are taken from piston nominal diameter. The circumferential traces shown in Figure 6.14(b) are measured 30.1752 mm down from skirt top 157 (bottom of oil ring groove). The ovalities are calculated from the circumferential trace. The vertical trace is combined with the piston ovality measurements to create the 3D piston skirt profile in CASE. For structural and thermal analysis, however, a complete 3D volume Computer-Aided Design (CAD) model is required. This is generated in NX and used in the meshing software Altair HyperMesh to generate linear tetrahedral mesh (see Figure 6.15). The mesh files are exported to CASE. The surface properties are obtained via a detailed analysis of the 3D optical profiler images using Omnisurf3D software (discussed in previous sections). 6.5.2 Block 2 (Running model using CASE-PISTON) Based on these inputs, the Cylinder-kit Assembly System for Engines (CASE) is used to model the cylinder kit components. CASE has two modules: CASE-RING and CASE-PISTON. In the current study, the model is set up in the CASE-PISTON module. The slider-crank mechanism is used to model the piston’s axial motion. The piston’s lateral and rotational motions are solved numerically when radial clearance exists between the piston and the cylinder bore under the analysis of piston elastohydrodynamics. CASE includes the lubrication model using Reynolds’ equation [41], asperity contact model using Greenwood-Tripp equation [43]. The competing forces acting on the piston such as combustion gas force, weight, and inertia of piston, wrist pin, connecting rod; crankshaft reaction force; hydrodynamic, contact, shear, frictional forces, etc. are considered in CASE-PISTON. Note that, for lubrication, CASE assumes a fully flooded condition on the piston skirt. 158 Figure 6.15 Piston mesh in CASE Table 6.8 Mesh details Mesh Size (mm) No. of Elements 5 14259 4.5 18072 Piston thermal deformation is not considered due to a lack of sufficient data to model an accurate temperature profile on the piston. The cylinder bore deformation is generally supplied as input data either from experimental measurements or from finite element analysis of the cylinder block. Since adequate data on bore deformation is not available, the piston thermal deformation is also neglected for consistency. Lubrication-induced skirt deformation is included. Guyan reduction [144] is used to obtain the compliance matrix for the piston skirt. The friction coefficient between the piston and cylinder is assumed to be 0.05 for the worn stock piston and 0.01 for the worn coated piston. Currently, no experimental friction coefficient data for the coated piston is available. Based on the range of variation of the surface properties listed in Table 6.4 and Table 6.5, it is assumed that the friction coefficient for the coated piston is one-fifth of the stock piston. This coefficient will be verified 159 experimentally in the future. The three-dimensional analysis in CASE-PISTON is conducted for both worn coated and worn stock pistons. The physical difference between the stock and coated pistons comes from the difference in the coating material composition and the way the coating is applied on the piston skirt. The APC is applied in thicknesses from 50µm to over 65µm on the diameter, so 25-33 µm per side going into the engine. On the other hand, the stock coating thickness is about 13 µm per side. Surface properties inputs such as - summit radius, summit density, and surface roughness are different for worn stock and worn coated models as observed in previous sections. The piston skirt profiles are also different, see Figure 6.14. For both worn coated and worn stock pistons, linear tetrahedral mesh files from HyperMesh are used. Although they are by nature stiff elements, the current version of CASE only allows this specific mesh type. In the future, this feature will be modified to incorporate other mesh element types. The study is conducted for two mesh sizes. Details are in Table 6.8. Both mesh sizes resulted in identical results. Therefore, the results from the 5mm mesh study are presented in the following sections. The three-dimensional analysis incorporates the implementation of variable surface properties such as summit density, summit radius, and surface roughness in the Greenwood-Tripp model. Figure 6.16 Surfaces showing cylinder-coated piston skirt and cylinder-stock piston skirt 160 Figure 6.16 shows the cylinder- piston skirt surfaces considered in Greenwood-Tripp model for both coated and stock piston skirts. Greenwood-Tripp model provides an expression for the nominal pressure carried out by surface asperities. Contact pressure is proportional to the square of the product of summit radius, summit density and roughness. 16√2π σ h pa (h) = (σs βη) 2 E ∗ √ s F5 ( ) (6.1) 15 β 2 σs Where, σs is combined surface roughness, β is average summit radius, η is surface density of asperity peaks and E ∗ is the effective modulus expressed in terms of the modulus and viscosity of each surface.h is the effective clearance between the skirt and liner. It is assumed that this effective clearance is equal to the oil film thickness which is calculated from the nominal piston to bore clearance, skirt deformation, contribution to oil film thickness due to piston ovality and skirt profile height. h is measured as the gap between the reference planes of each surface. Some surface asperities are above this reference plane while some are below. The interaction of the asperities above the reference plane is what results in contact friction. In short, while the reference planes are separated by some distance, the surface asperities still touch. There is always oil film thickness while asperity contacts are still engaging. It does not require that the oil film thickness be zero for contact to engage. This is what is known as the mixed-lubrication regime. Contact engages when effective clearance between the piston-cylinder surfaces (h) is less or equal to 4 σs (roughness)- h ≤ 4 σs contact engages h > 4 σs hydrodynamic lubrication (6.2) Figure 17(a) and Figure 17(b) shows the asperity density for worn stock and worn coated piston 161 skirt, respectively. A larger area with higher values around the centerline on the worn coated piston skirt is observed which would help the surface to retain more oil and reduce friction. Figure 18(a) and Figure 18(b) shows the asperity radius for worn stock and worn coated piston skirt, respectively. It is observed that at the center location, the radius is higher for the worn coated skirt than the worn stock. A larger area with higher values around the centerline of the worn coated piston skirt is seen. At the side locations, both have a radius of 2.74 µm to 25 µm. On the worn coated skirt, the radius is higher at the top and lower at the bottom. On the worn stock skirt, the radius is lower at the bottom and higher at the top. Figure 6.17 Piston skirt asperity density (a) worn stock (b) worn coated. 162 Figure 6.18 Piston skirt asperity radius (a) worn stock (b) worn coated. Figure 6.19 Piston skirt surface roughness (a) worn stock (b) worn coated. Figure 19(a) and Figure 19(b) shows the surface roughness for worn stock and worn coated piston skirt, respectively. It is observed that roughness is significantly higher on the worn coated piston skirt than on the worn stock piston skirt. This means higher oil retention and reduced friction for the worn coated piston skirt. 163 6.6 CASE Simulation Results Results comparing the worn stock and worn coated piston skirt performance are discussed in this section. All the plots presented in this section are averaged across the surface. Throughout the cycle, the piston is repeatedly forced up against the cylinder liner. The force from the piston is supported by two reaction forces. There is a hydrodynamic pressure force from the lubricant and an asperity contact force, which occurs when the asperities of two surfaces start to directly interact. The hydrodynamic pressure profile is determined via Reynold’s equation [41] which does not consider the effect of surface properties. The asperity contact pressure is determined via the Greenwood -Tripp model [43] that accounts for the effect of surface properties such as summit radius, summit density, and surface roughness. As the piston travels laterally within the cylinder bore, the oil film thickness can become very small, thus allowing for possible solid-to-solid contact. The magnitude of the pressure arising by this contact is calculated via the Greenwood-Tripp model. Integrating the asperity contact pressure over the contact surface yields the total load carried out by the asperities. Figure 6.20 shows the averaged surface-to-surface contact forces comparison between the coated and stock models. Given the surface measurements, a worn coated piston should have higher contact forces. It is observed that the contact forces are higher for the coated model at both major and minor thrust sides. 164 Figure 6.20 Contact force for coated vs stock piston skirt. Figure 6.21 Hydrodynamic force comparison between the worn coated & worn stock piston skirts. 165 The hydrodynamic pressure is calculated using Reynold’s equation. Once the hydrodynamic pressure profile across the cycle is found, it can be integrated across the surface to find the total hydrodynamic pressure force. Figure 6.21 presents the hydrodynamic forces comparison between the coated and stock models. It is observed that the hydrodynamic force is very small for both models compared to contact force. Note that the hydrodynamic force is calculated using classical Reynold’s equation which does not consider the effect of surface properties. These properties are the key factor distinguishing the coated piston skirt from the stock skirt. Figure 6.22 Total force comparison between the worn coated & worn stock piston skirts. The total force, which is the sum of the hydrodynamic force and contact force, shows the dominant effect of contact force, see Figure 6.22. Both curves overlap each other at most crank angles. The oil film thickness on the piston skirt depends on the nominal piston to cylinder bore clearance, cylinder bore deformation, skirt deformation, contribution to oil film thickness due to piston ovality, skirt profile height, and skirt eccentricity at the top, bottom, and along the wrist pin [145]. 166 In this study, the cylinder bore deformation and corresponding piston thermal deformation are neglected. Figure 23(a) and Figure 23(b) shows the minimum oil film thickness averaged across the skirt surface comparison between the worn coated & worn stock piston skirt at major thrust side and minor thrust side, respectively. At both sides, the coated piston skirt surface consistently maintains higher minimum oil film thickness across the cycle. Figure 6.23 Minimum skirt oil film thickness comparison between the worn coated & worn stock piston skirts (a) At major thrust side (b) At minor thrust side. As the piston slides relative to the cylinder bore, it experiences friction resisting its motion. The total friction is due to both the hydrodynamic friction which results from the shearing of a viscous lubricant, as well as the contact frictional force which results from the interacting surface asperities. 167 Figure 6.24 Hydrodynamic shear force comparison between the worn coated & worn stock piston skirts. Figure 6.24 shows the hydrodynamic shear force comparison between the worn coated & worn stock piston skirts. Hydrodynamic shear stress (τ) is given by 𝜇𝜈𝑝 ℎ 𝜕𝑃ℎ 𝜏= − +2 (6.3) ℎ 𝜕𝑦𝑠 Here, 𝜇 is dynamic viscosity, 𝜈𝑝 is piston axial velocity, ℎ is oil film thickness, 𝑃ℎ is hydrodynamic pressure and 𝑦𝑠 is skirt axial location. The second term in equation 6.3 is negligible relative to the first term. So, the hydrodynamic shear stress is inversely proportional to the oil film thickness. It is observed from Figure 6.24 that, the hydrodynamic shear force is consistently lower across the cycle for the coated piston skirt. This also corroborates the findings from Figure 6.23 that, the coated piston skirt always maintains a higher minimum oil film thickness. 168 Figure 6.25 Contact friction force comparison between the worn coated & worn stock piston skirts. Figure 6.26 Total friction force comparison between the worn coated & worn stock piston skirts. It is observed that the friction forces are higher on the stock piston than the coated one only during expansion stroke, see Figure 6.25. The coated friction force is higher rest of the cycle. Note that, 169 friction coefficients used for the coated and stock piston skirts are 0.01 and 0.05, respectively. Figure 6.26 shows the total friction force which is the sum of hydrodynamic shear force and contact friction force. This also follows the same trend as Figure 6.25. The friction power loss comes from the net friction work on the skirt due to excessive contact with the cylinder bore. This is undesirable, as it can result in the net power loss of the engine and increased wear of the skirt and cylinder bore surfaces. The average friction power loss and net friction work on the skirt per cycle are calculated and listed in Table 6.9 for both conditions. The worn coated piston skirt resulted in a slightly higher friction power loss compared to the worn stock piston. This is also demonstrated in Figure 6.27 that presents the total friction power loss. It is observed that the coated piston skirt experiences a lower friction power loss during expansion stroke and higher during rest of the cycle. Note that friction coefficients used for the coated and stock piston skirts are 0.01 and 0.05, respectively. Table 6.9 Performance comparison between worn coated and worn stock conditions Performance Parameter Worn Worn Worn Worn Worn coated coated stock coated coated 25% original original 50% reduced skirt skirt reduced skirt skirt Friction coefficient 0.01 0.05 0.05 0.01 0.01 Average friction power loss (W) 2169.6 1820.9 10716 255.3 596.27 Net friction work on skirt per 128.91 97.61 637.18 13.664 33.454 cycle (J) 170 Figure 6.27 Friction power loss comparison between the worn coated & worn stock piston skirt Figure 6.28 Friction power loss using friction coefficient 0.05 for both stock and coated piston model Figure 6.28 presents an extreme case where the coated piston skirt friction coefficient is 0.05, same 171 as the friction coefficient of the stock piston skirt. The friction power loss at this condition is predicted to be very high, see Table 6.9. To explore the effect of coated piston skirt profile variation, the coated skirt profile is reduced by 25% and 50%. Two simulations are conducted keeping everything else the same as the original coated piston model. Friction coefficient 0.01 for coated pistons and 0.05 for stock piston model are used. Figure 6.29 shows significant reduction in piston skirt profile as the coated piston skirt profile is reduced. The corresponding friction power loss and net friction work are listed in Table 6.9. Figure 6.29 Friction Power Loss for different projected piston skirt profiles 172 Figure 6.30 Piston secondary motion comparison between the worn coated & worn stock piston skirts (a)piston tilt (b) piston eccentricity (translation) perpendicular to the wrist pin The transverse components of the combustion gas force and inertia arising due to connecting rod orientation result in a translation motion perpendicular to the piston pin axis and rotation motion about the piston pin. The translation motion is also called lateral motion and the rotational motion is called piston tilt. Piston tilt is calculated based on three parameters- piston top eccentricity, or the transverse distance from the center point on top of the piston to the center of the bore axis, piston bottom eccentricity, or the transverse distance from the center point on the bottom of the piston to the center of the bore axis and the piston height. These two motions constitute the piston secondary motion. The hydrodynamic and contact forces and moments, as well as the shear and contact friction forces and moments, are functions of the piston secondary motion [145]. Figure 6.30(a) and Figure 6.30(b) show that the coated piston experiences lower piston tilting and eccentricity at the wrist pin level. 6.7 Summary and Conclusions In this study, the unique surface characteristics of APC coated, and stock piston skirts have been investigated and those properties varying across the surface have been used to develop a numerical 173 model in CASE-PISTON. This research will advance the understanding of the detailed physics of cylinder-kit performance and complex interactions between lubrication and friction of the coated pistons in comparison to the stock pistons. This will enable to quantify the influence of applying the abradable powder coating to a piston assembly. The surface images of the worn coated and stock piston skirts have been generated using a 3D optical profilometer. Image analysis has been conducted using Omnisurf3D. The surface properties show that as the surface wears the smaller structures are crushed leaving a higher summit radius for both coated and stock piston skirts. The coated piston skirt surface has significantly higher surface roughness than the stock piston skirt which makes it a better candidate for oil retention and lubrication. After wearing, the surface is left with numerous microscopic ripples that result in higher summit density. The valley structure is believed to play a significant role in the coated surface in terms of oil retention and lubrication. It is observed that the coated surface has deeper valleys. The variable surface property distribution shows a larger area with higher values of summit radius and summit density around the centerline on the worn coated piston skirt. This would help the surface to retain more oil and reduce friction. CASE results are presented to demonstrate the capability of the modeling methodology to capture the unique variable surface properties of the coated and stock pistons. The results show that the cylinder-piston skirt contact force is higher for the coated piston skirt than the stock ones based on the Greenwood-Tripp model. The cylinder-piston skirt contact force is found to be higher for the coated piston skirt than the stock ones using the Greenwood-Tripp model. There is always oil film thickness while asperity contacts are still engaging. It does not require that the oil film thickness be zero for contact to engage. The total force which is the sum of the hydrodynamic force and contact force shows dominant effect of contact force. The minimum oil film thickness of the coated 174 piston skirt is always higher which results in consistently lower hydrodynamic shear force for the coated piston skirt. The friction forces and friction power loss show that the coated piston skirt experiences higher friction power loss even when friction coefficient used for the coated and stock piston skirts are 0.01 and 0.05, respectively. These coefficients must be verified. Using the same friction coefficients for both results in even higher friction power loss for the coated piston skirt. Reducing the coated piston skirt profile shows significant improvement in terms of friction power loss. The coated piston experiences lower piston tilting. The hydrodynamic force is calculated using Classical Reynold’s equation that does not consider the effect of surface properties which is the key factor distinguishing the coated piston skirt from the stock skirt. The classical Reynold’s equation appears to be inadequate to address the effect of surface microstructure on hydrodynamic pressure development because the surface roughness is not being considered. In future works, the modified Reynold’s equation developed by Patir and Cheng [42] will be implemented to include the surface roughness effect in the effort of modeling piston assembly lubrication. This equation will account for the effect of surface roughness with the addition of flow factors. Another flow factor to be added in the modified Reynold’s equation will account for the valleys that retain oil. A study will be conducted to explore how the change in piston tilt due to coating is affecting the performance of rings. The results presented here are obtained by running the model for 1 cycle. In the future, the simulation will be conducted for 1000s of hours to explore the effect on wear and wear coefficients. Simulations will be conducted for a range of operating conditions. 175 Chapter 7 Concluding Remarks and Recommendations for Future Work 7.1 Coated Piston Performance Assessment: Concluding Remarks and Recommendations for Future Work The study focuses on surface properties characterization of Line2Line novel abradable powder coated piston skirts and stock piston skirts of a Cummins 2.8 L turbo engine. The surface images of the worn coated and stock piston skirts have been generated using 3D optical profilometer and image analysis has been conducted using Omnisurf3D. Later these properties along with other cylinder-kit properties are used to build a model in CASE. During surface properties measurement, only five measurements on each surface are taken and the rest are interpolated or extrapolated during the implementation of variable surface properties in the Greenwood-Tripp Model. In future more locations at each surface should be measured on each surface for better representation of the surfaces. In the current model it is assumed that the piston skirt is fully flooded. There is always oil film thickness while asperity contacts are still engaging. It does not require that the oil film thickness be zero for contact to engage. Since for the performance assessment of the L2L coating and stock coating we are concerned about which coated surface can retain more oil, starved lubrication implementation should be part of future CASE PISTON program development. In the model the friction forces and friction power loss are calculated using friction coefficient of 0.01 for the coated and 0.05 for the stock piston skirts. The surface properties of the coated surface 176 are observed to be 1.5-9 times higher at different locations according to the surface measurement. Based on this the coated piston skirt friction coefficient is considered to be 1/5 of that of the stock piston skirt. Hypothetically considering same friction coefficient for both models demonstrated that, the APC coated piston skirt experienced higher friction power loss. Therefore, these coefficients must be verified. In the future, the coefficient between coated piston skirt and cylinder will be verified through a combination of tribometer testing and numerical model calibration. A separate model on sub-micron level can be developed replicating a cylinder liner and a piston skirt sliding along it. The friction coefficient can be estimated by applying all the forces and lubrication properties of the surfaces. The piston coating used in the current study is part of the first iteration of the process. The coating appeared to be too thick for the skirt-liner clearance. Reducing the coated piston skirt profile in the model shows significant improvement in terms of friction power loss. A study should be conducted using CASE and HEEDS to find the coated skirt profile that would optimize both friction power loss and piston tilt. In the current version of CASE-PISTON, classical Reynold’s equation is used to model the hydrodynamic pressure between the cylinder liner-piston skirt surface. The classical Reynold’s equation appears to be inadequate to address the effect of surface microstructure on hydrodynamic pressure development because the surface roughness is not being considered. The modified Reynold’s model developed by Patir and Cheng [42,146,147] can be adopted to include surface roughness effect in the effort of modeling piston assembly lubrication. The governing equation for the model is a modified version of the classical Reynold’s equation but it includes shear and pressure flow factors [148,149], which are functions of the ratio between surface roughness and nominal clearance between two surfaces. These flow factors are used for the consideration of the 177 effect of surface characteristics on the hydrodynamic flow and pressure generation. Basically, they are dependent on the surface microstructure and the gap between two relatively moving surfaces. The Patir and Cheng model is already implemented in CASE-RING program, and currently being implemented in the PISTON program of CASE. Another factor to be added in the modified Reynold’s equation will account for the valleys which retain oil. A study should be conducted to explore how the change in piston secondary motion due to coating is affecting the performance of rings. A slight change in the piston tilt angle affects the hydrodynamic lubrication of the piston and rings. Piston tilt angle can cause rings to switch from contact forces to hydrodynamic forces, thus reducing friction. Oil film thickness left on the cylinder liner by the piston/rings could be different for the two pistons. The current version of CASE considers uniform ring twist across the circumference. But the 3D ring FEA contact model has shown that the twist angle varies across the circumference. In future, non-uniform ring twist should be incorporated in CASE. 7.2 3D Multiphysics Methodology: Concluding Remarks and Recommendations for Future Work The current study is the first modeling effort toward developing a three-dimensional, multi-physics methodology to investigate liquid oil-oil vapor and combustion gas transport in the cylinder-kit assembly. The results are specific to the cylinder-kit geometry and engine operating conditions. First, a 50 mm small bore engine is considered to implement the methodology. The methodology is applied for a non-twisted and a negatively twisted second ring cylinder-kit configurations. The key findings of the model show oil and gas mass distribution in all parts of the domain and their impact on blow-by, reverse blow-by and oil consumption. The results are within the acceptable range for typical engines. A detailed mesh sensitivity study indicated the need for a very fine mesh 178 size for capturing the oil-air interface as well as certain areas where the ring-groove or ring-liner clearance is considerably low. Table 7.1 Simulation time and grid details of the cases Case Cell size Total no. of No. of # 𝐜𝐞𝐥𝐥𝐬 Total CPU Total runtime no. (µm) cells Cores # 𝐜𝐨𝐫𝐞𝐬 hours 1 500 25,238 50 505 750 15 hrs 2 250 171,435 100 1715 74,800 1080 hrs ≈ 1.5 months 3 125 1,266,257 250 5066 443,000 1872 hrs ≈ 2.6 months 4 100 2,403,538 250 9615 522,000 2188 hrs≈3 months 5 75 4196137 300 13988 760,000 2640 hrs≈ 4 months Table 7.1 shows the simulation time for each grid size. It is observed that the real bottle neck of the methodology is the simulation time. For a grid size of 75 µm (see Case no. 5), no difference in simulation speed is observed using 300 cores or 500 cores, This indicates that scalability is the main issue. In other words, adding more core after a certain point does not help. In High performance computing (HPC) scalability or scaling is widely used to indicate the ability of hardware and software to deliver greater computational power when the number of resources is increased. This depends on the exact hardware configuration and the details of the simulation. There were some additional blockades like system scheduling limitation that added to the already “runtime constrained” simulation. At MSU HPC, each user can use a maximum of 500 cores. The scheduler also attempts to ensure fair resource utilization of all HPCC users by adjusting the initial priorities of the users who have recently used HPCC resources. Due to the policy, if users had jobs running with many resources recently, their current pending jobs might wait longer than before. The comparison between Intel MPI vs. OPENMPI showed that Intel MPI ensured faster 179 simulation. There are also file system quota issues inherent to the system. MSU HPC is still working to fix it. In addition to these, restarting and mapping files are corrupted often during the runtime, which also adds to the total run time. Figure 7.1 The three rings oil film thickness prediction by CASE Due to the long run time several simulation features are suppressed. They are- • All aspects of cylinder-kit motion are not directly incorporated in the CONVERGE model. The piston and ring motion, ring twists and resultant ring flutter or collapse (if any) are considered in the ring FEA contact model and consequently, the geometry of the cylinder- kit is dynamically varied throughout the cycle in the CONVERGE CFD model. But the effect of inertia is not directly incorporated in the CFD model as a piston or ring motion boundary condition. • It is observed that there are regions at the ring-groove, ring -liner interfaces where the clearance is as low as 2.5 microns. Figure 7.1 shows the oil film thickness at the ring-liner interface, which indicates such narrow clearance. Fixed embedding on top ring and second 180 rings to ensure finer mesh close to walls is required in those regions. The embedded cell size depends on the base size based on equation (7.1) 𝐵𝑎𝑠𝑒 𝑠𝑖𝑧𝑒 Embedded cell size =2𝑒𝑚𝑏𝑒𝑑 𝑠𝑐𝑎𝑙𝑒 (7.1) A 75 µm base grid with embed scale of 5 results in an embedded cell of 2.34 µm which would increase the wall clock time significantly. Hence, this feature is not incorporated. • Adaptive meshing is also not implemented for the same tradeoff among adaptive mesh refinement (AMR) sub-grid scale embedding level, sub-grid criterion (0.1-1% of the characteristic value) and wall clock time. The experimental validation of the methodology is conducted on a 137.02 mm bore Cummins Acadia engine. The Table 7.2 shows a mesh projection. The nominal piston-bore clearance is 45 µm. Table 7.2 shows that a grid size of 45 µm would result in a 124.81 M cells which would require years to complete. So, based on previous simulation runtime the 150 µm mesh size was selected that took about four months to complete the simulation. Table 7.2 Mesh size and no of cells projection for Acadia engine Case no. Cell size (µm) Total no. of cells 1 200 1.42M 2 150 3.9M 3 125 5.82M 4 100 11.36M 5 75 27M 6 45 124.81M To better understand the real bottleneck of the simulation a strong scaling comparison is conducted between MSU HPCC and Argonne THETA. THETA uses “Intel Xeon Phi x200, codename 181 Knights Landing (KNL)”, a second-generation many integrated core (MIC) architecture product from Intel and MSU HPC uses “Intel Xeon @ 2.4 GHz+ and AMD EPYC @2.5 GHz+”. A single core of one of these Xeons can be 5-6x faster than a single KNL core. So, we need to use many more KNLs to exceed the performance of a collection of Xeons. The one-hour runtime strong scaling experiment shows that adding number of CPUs is not speeding up the simulation in either of the computing facility. MSU HPC is found to be faster. Firstly, the Xeon Gold cores are faster than Phi cores used at Argonne. Secondly, MSU HPC is optimized for throughput so that it can speed up smaller jobs. The latency comparison shows that at MSU the lowest is 3 and highest is 12 microseconds. At Argonne, the lowest is 8 microseconds. So, it appears that, Argonne’s network is slightly slower. Note that, latency varies based on job load and trials. The real bottle neck could be memory bound, network bound or both. Memory bandwidth is the rate at which CPUs can read data from RAM. Again, the network communication between nodes could be affecting the runtime as well. Here the ability of the problem/code to scale is key. The current problem configuration is not obtaining load balancing well in CONVERGE. No matter how many CPUs are added, after a certain point the speed is not increasing. Either the problem configuration needs to be reformulated or simplified. Or a different solver/algorithm needs to be used. Potential approaches to solve the “simulation time” issue are listed below- • Using a solver that supports GPUs. The new generation GPUs scale well across nodes and have higher memory bandwidth. But CONVERGE does not support GPUs. • Another solution can be problem specific heuristic understanding and development of a 182 custom solver. No commercial solvers might be able to solve this specific complex problem completely. • Another approach can be a “Hierarchical Solution”. First the problem will be solved for a coarse mesh size. Then using the macro level information from the coarser mesh study, finer mesh studies can be conducted at certain range of crank angles. Some model specific near-term improvements that should be considered are- • Considering piston and bore deformation effect in both ring FEA contact model and LINCC program. Currently, we assume that both piston and bore are undeformed. • Currently the LINCC program can generate dynamically varying geometries in STL format by combining the rings, piston, and liner. But the geometry contains a significant number of overlapping or intersecting triangles that need to be cleaned before using in the CONVERGE CFD solver. The geometry generation algorithm should be modified to account for these error triangles to ensure the generation of clean geometries to be fed directly to the CONVERGE CFD solver. • Evaluation of inertia effect by directly incorporating piston/ring motion in the model if the computational time issue can be solved. • In the current methodology, the oil is initialized as a mass fraction. A better approach would be taking the oil film estimate from CASE and initializing the oil film profile across the domain. A separate oil filling ratio variation and oil re initialization at the piston skirt before the expansion stroke indicated that the prescribed oil is entirely problem specific. Further studies are required to ensure optimum oil supply in the domain across the cycle. • Model calibration based on experimental validation- 183 The experimental validation showed that the model predicts 3-7 times higher values for land pressures, blowby and typical engine oil consumption. There is substantial amount of variability in the experimental results. The experimental land pressures are measured only at a certain location that might not be representative of the other regions. Again, the modeling methodology has used a second ring (Part 3685672) that has a different ring end gap than the second ring (Part 2869635) used in the experiment. Part 3685672 was used for model since the ring Part 2869635 was not available in the market and Cummins could not provide with the ring free shape profile. Figure 7.2 shows how the twist angle variation be totally different for these two rings. This means, there would be significant differences between the geometries created using the two different rings. Figure 7.2 Difference of the twist angle circumferential variation between the two different ring parts at a selected CA Ring twist analysis via the ring FEA contact model is highly critical in the result prediction of the methodology. In future, experimental validation and subsequent model calibration should be carried out using the exact same geometry parts. • A complete evaporation model sensitivity study with appropriate hydrocarbon species- 184 • The current methodology assumes n-heptane to be the lightest constituent, since the exact composition is not known. The exact oil composition should be incorporated in collaboration within oil company. Also, a complete evaporation model sensitivity study needs to be conducted. • Incorporating a feedback loop where the results of the 3D model will be used as inputs for CASE model to reevaluate the piston/ring motion and other features should be considered. In order to improve the accuracy of the model, the ring dynamics needs to be coupled with the 3D gas dynamic model in an iterative manner. This means not only the ring motion will be affected by the pressure distribution; the local gas pressures will also be influenced by the ring motion. Computational power will change over the next decade, giving a great advantage to those that employ this methodology. In the presence of a large computing platform the long-term goal is to use the 3D multiphysics modeling methodology to optimize designs for low friction, low oil consumption and long life. If successful, the new methodology could completely disrupt the design of cylinder-kits. Ideally, we would like to have this methodology produce representative designs within hours or a week. This would promote the option of conducting major optimizations based on the performance of proposed designs. Specifically, optimization of a particular design will require implementation of very efficient “supervisory optimization” codes. One potential example is the implementation of HEEDS, which is used in chapter 4 with CASE. Figure 7.3 presents a. flowchart showing cylinder-kit design optimization using 3D multiphysics methodology. 185 Figure 7.3 A flowchart showing cylinder-kit design optimization using 3D multiphysics methodology. The key challenge in bringing this methodology to the market is that the computational capability will not evolve rapidly enough to make 3D modeling an everyday tool. To tackle this, using a hierarchical solution with the combination of a high-level coarse mesh model and multiple fine mesh models at critical crank angles/regions of the domain across the cycle can be an option. Currently, a complete 3D model takes about 4- 5 months from start to end whereas the CASE model only takes few minutes to complete. Shorter simulation time means faster iteration. So, the lessons learned from the 3d methodology can be used to improve current CASE simulation results. Flow coefficients at the ring-liner interface, ring-groove interface or ring end gap can be calculated from the 3D model to better estimate the blowby, ring motion using the CASE model in a short 186 time. This methodology addresses development of a design tool that can influence the design of hundreds of millions of cylinder-kit assembly components manufactured annually.The manufacturers of internal combustion engines and components can benefit from this effort. These companies provide pistons, piston rings, and coatings for cylinder-kit components. Worldwide, approximately 200 million new internal combustion engines are manufactured annually, each with one or more cylinder-kit. With replacement pistons and cylinder-kits for applications other than IC engines, there may be 1 trillion cylinder-kit assemblies manufactured annually. The potential societal value of the methodology includes more effective, long-lived aftertreatment systems and fuel savings. The new design system will promote opportunities for new products, such as an advanced ring configurations and cost-effective micron-coatings of select cylinder-kit parts. Thus, the methodology will elucidate new physics, promote new cylinder-kit designs and has the potential to influence manufacturing of these systems on a grand scale. 187 APPENDIX 188 Figure A-1. (a) 3D contour plot at side location 1 (A1) of the coated post run worn piston skirt (b) Zoomed in view of the contour plot (c) Location A1 labeled on the coated post run worn piston skirt 189 Figure A-2. (a) 3D contour plot at center location 1 (B1) of the coated post run worn piston skirt (b) Zoomed in view of the contour plot (c) Location B1 labeled on the coated post run worn piston skirt Figure A-3. (a) 3D contour plot at side location 1 (C1) of the coated post run worn piston skirt (b) Zoomed in view of the contour plot (c) Location C1 labeled on the coated post run worn piston skirt 190 Figure A-4. (a) 3D contour plot at center location 1 (D1) of the coated post run worn piston skirt (b) Zoomed in view of the contour plot (c) Location D1 labeled on the coated post run worn piston skirt 191 Figure A-5.GT map of the Cummins R2.8L turbocharged four-cylinder diesel engine 192 BIBLIOGRAPHY 193 BIBLIOGRAPHY 1. Atis, C., Chowdhury, S.S., Ayele, Y., Stuecken, T., Schock, H., and Voice, A.K., “Ultra- Lean and High EGR Operation of Dual Mode, Turbulent Jet Ignition (DM-TJI) Engine with Active Pre-chamber Scavenging,” SAE Technical Paper, ISBN 0148-7191, 2020. 2. Richardson, D.E., “Review of power cylinder friction for diesel engines,” Proc. 1999 Spring Tech. Conf. ASME Intern. Combust. Engine Div. (Volume 3), Columbus, Indiana, USA, 24- 28 April 1999 122(October):71, 1999, doi:10.1115/1.1290592. 3. Ingwersen, J., Pederson, P.S., Nielsen, T., Larsen, E., and Fenger, J., Effects on Pollution of a Reduction or Removal of Lead Addition to Engine Fuel, 1978. 4. Trier, C.J., Coates, D.S., Fussey, D.E., and Rhead, M.M., “The contribution of lubricating oil to diesel exhaust emissions,” Sci. Total Environ. 93:167–174, 1990. 5. Sakurai, H., Tobias, H.J., Park, K., Zarling, D., Docherty, K.S., Kittelson, D.B., McMurry, P.H., and Ziemann, P.J., “On-line measurements of diesel nanoparticle composition and volatility,” Atmos. Environ. 37(9–10):1199–1210, 2003. 6. Ruddy, B.L., Dowson, D., Economou, P.N., and Baker, A.J.S., “Piston-ring lubrication. Part III. The influence of ring dynamics and ring twist,” Energy Conserv. through Fluid Film Lubr. Technol. Front. Res. Des. 191–215, 1979. 7. Tian, T., Rabute, R., Wong, V.W., and Heywood, J.B., “Effects of piston-ring dynamics on ring/groove wear and oil consumption in a diesel engine,” SAE Tech. Pap. (412), 1997, doi:10.4271/970835. 8. Tian, T. and Wong, V.W., “Modeling the lubrication, dynamics, and effects of piston dynamic tilt of twin-land oil control rings in internal combustion engines,” J. Eng. Gas Turbines Power 122(1):119–129, 2000. 9. Tian, T., “Dynamic behaviours of piston rings and their practical impact. Part 1: ring flutter and ring collapse and their effects on gas flow and oil transport,” Proc. Inst. Mech. Eng. Part J J. Eng. Tribol. 216(4):209–228, 2002. 10. Thirouard, B. and Tian, T., “Oil transport in the piston ring pack (Part I): identification and characterization of the main oil transport routes and mechanisms,” SAE Technical Paper, ISBN 0148-7191, 2003. 11. Thirouard, B. and Tian, T., “Oil transport in the piston ring pack (Part II): zone analysis and macro oil transport model,” SAE Technical Paper, ISBN 0148-7191, 2003. 12. Thirouard, B., Tian, T., and Hart, D.P., “Investigation of oil transport mechanisms in the piston ring pack of a single cylinder diesel engine, using two dimensional laser induced fluorescence,” SAE Trans. 2007–2015, 1998. 194 13. Yilmaz, E., “Sources and characteristics of oil consumption in a spark-ignition engine.,” DSpace@MIT. Thesis (Ph.D), 2003. 14. Baelden, C. and Tian, T., “A dual grid curved beam finite element model of piston rings for improved contact capabilities,” SAE Int. J. Engines 7(1):156–171, 2014. 15. Liu, Y., Li, Y., and Tian, T., “Development and Application of Ring-Pack Model Integrating Global and Local Processes. Part 2: Ring-Liner Lubrication,” SAE Int. J. Engines 10(4), 2017, doi:10.4271/2017-01-1047. 16. Bhouri, M.A., Curved beam based model for piston-ring designs in internal combustion engines, 2017. 17. Wang, Y., Air flow effects in the piston ring pack and their implications on oil transport, 2012. 18. Fang, T., Computations and modeling of oil transport between piston lands and liner in internal combustion engines, 2014. 19. Przesmitzki, S.S.V., Characterization of oil transport in the power cylinder of internal combustion engines during steady state and transient operation, 2008. 20. Felter, C.L., “Lubrication of Piston Rings in Large 2 – and 4 – stroke Diesel Engines,” 2007. 21. Shahmohamadi, H., Rahmani, R., Rahnejat, H., King, P., and Garner, C., “Cavitating flow in engine piston ring-cylinder liner conjunction,” ASME 2013 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers Digital Collection, 2013. 22. Shahmohamadi, H., Mohammadpour, M., Rahmani, R., Rahnejat, H., Garner, C.P., and Howell-Smith, S., “On the boundary conditions in multi-phase flow through the piston ring- cylinder liner conjunction,” Tribol. Int. 90:164–174, 2015. 23. Hronza, J. and Bell, D., “A lubrication and oil transport model for piston rings using a Navier-Stokes equation with surface tension,” SAE Tech. Pap. 2007(724), 2007, doi:10.4271/2007-01-1053. 24. Puthiya Veettil, M. and Shi, F., “CFD analysis of oil/gas flow in piston ring-pack,” SAE 2011 World Congr. Exhib., 2011, doi:10.4271/2011-01-1406. 25. Oliva, A., Held, S., Herdt, A., and Wachtmeister, G., “Numerical Simulation of the Gas Flow Through the Piston Ring Pack of an Internal Combustion Engine,” SAE Technical Paper, ISBN 0148-7191, 2015. 26. Oliva, A. and Held, S., “Numerical multiphase simulation and validation of the flow in the piston ring pack of an internal combustion engine,” Tribol. Int. 101:98–109, 2016, doi:10.1016/j.triboint.2016.04.003. 195 27. Lyubarskyy, P. and Bartel, D., “2D CFD-model of the piston assembly in a diesel engine for the analysis of piston ring dynamics, mass transport and friction,” Tribol. Int. 104(September):352–368, 2016, doi:10.1016/j.triboint.2016.09.017. 28. Dr. P. Urzua, A. Lehnen, Dr. R. Otte, Dr. M. Prouvier, A.W.– V.A.P.D.T. and Tian, Q. Zhang, Y.L.– M., “Detailed Fuel and Oil Transport Simulation for Understanding Low Speed Pre-Ignition – LSPI,” 13th International AVL Symposium on Propulsion Diagnostics, Kurhaus, Baden-Baden,Germany, 2018. 29. CASE 6.0 Theoretical Manual, Mid Michigan Res., 2020. 30. Brombolich, L.J., “Structural Mechanics of Piston Rings,” Compu-Tec Eng. Intern. Report, St. Louis, MO, 1993. 31. Ting, L.L., “A review of present information on piston ring tribology,” SAE Tech. Pap., 1985, doi:10.4271/852355. 32. Ting, L.L. and Mayer Jr, J.E., “Piston ring lubrication and cylinder bore wear analysis, part I—theory,” 1974. 33. Ejakov, M.A., Schock, H.J., Brombolich, L.J., Carlstrom, C.M., and Williams, R.L., “Simulation Analysis of Inter-Ring Gas Pressure and Ring Dynamics and Their Effect on Blow-by,” ASME Pap., 1997. 34. Ejakov, M.A., Hascher, H.G., Schock, H.J., Brombolich, L.J., Trinker, F.H., and LoRusso, J.A., “Simulation analysis of ring pack behavior in a deactivated cylinder,” American Society of Mechanical Engineers, New York, NY (United States), 1996. 35. Ejakov, M.A., Schock, H.J., and Brombolich, L.J., “Modeling of ring twist for an IC engine,” SAE Trans. 2284–2298, 1998. 36. Chui, B.-K., Schock, H.J., Fedewa, A.M., Richardson, D.E., and Shaw, T., “Coupled Model of Partially Flooded Lubrication and Oil Vaporization in an Internal Combustion Engine,” ASME 2005 Internal Combustion Engine Division Spring Technical Conference, American Society of Mechanical Engineers Digital Collection: 447–456, 2005. 37. Panayi, A., Schock, H., Chui, B.-K., and Ejakov, M., “Parameterization and FEA approach for the assessment of piston characteristics,” SAE Technical Paper, ISBN 0148-7191, 2006. 38. Panayi, A.P. and Schock, H.J., “Avenues for predicting piston wear: employing 2D and 3D numerical piston dynamics models,” SAE Int. J. Engines 1(1):713–722, 2009. 39. Cheng, C., Kharazmi, A., Schock, H., Wineland, R., and Brombolich, L., “Three- dimensional piston ring–cylinder bore contact modeling,” J. Eng. Gas Turbines Power 137(11), 2015. 40. Kharazmi, A., “Three Dimensional Analysis of the Gas Flow in Piston Ring Pack,” Michigan State University, ISBN 0355361442, 2017. 196 41. Reynolds, O., “IV. On the theory of lubrication and its application to Mr. Beauchamp tower’s experiments, including an experimental determination of the viscosity of olive oil,” Philos. Trans. R. Soc. London (177):157–234, 1886. 42. Patir, N. and Cheng, H.S., “An average flow model for determining effects of three- dimensional roughness on partial hydrodynamic lubrication,” 1978. 43. Greenwood, J.A. and Tripp, J.H., “The contact of two nominally flat rough surfaces,” Proc. Inst. Mech. Eng. 185(1):625–633, 1970. 44. Archard, J.F., “Contact and rubbing of flat surfaces,” J. Appl. Phys. 24(8):981–988, 1953, doi:10.1063/1.1721448. 45. Antoine, C., “Vapor Pressure: a new relationship between pressure and temperature,” Comptes Rendus Des Séances l’Académie Des Sci. (in French), 1888. 46. CONVERGE, C.F.D., “Software,(version 2.1. 0), convergent science,” Inc., Middleton, WI, 2013. 47. Manninen, M., Taivassalo, V., and Kallio, S., On the mixture model for multiphase flow, 1996. 48. Hirt, C.W. and Nichols, B.D., “Volume of fluid (VOF) method for the dynamics of free boundaries,” J. Comput. Phys. 39(1):201–225, 1981. 49. Converge, “CONVERGE Manual v2.4,” 1008, 2017. 50. Waclawczyk, T. and Koronowicz, T., “Modeling of the wave breaking with CICSAM and HRIC high resolution schemes,” Eur. Conf. Comput. Fluid Dyn. 1–19, 2006. 51. Fefferman, C.L., “Existence and smoothness of the Navier-Stokes equation,” Millenn. Prize Probl. 57:67, 2006. 52. Launder, B.E. and Sharma, B.I., “Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc,” Lett. Heat Mass Transf. 1(2):131–137, 1974. 53. Redlich, O. and Kwong, J.N.S., “On the thermodynamics of solutions. V. An equation of state. Fugacities of gaseous solutions.,” Chem. Rev. 44(1):233–244, 1949. 54. Issa, R.I., “Solution of the implicitly discretised fluid flow equations by operator-splitting,” J. Comput. Phys. 62(1):40–65, 1986. 55. Courant, R., Friedrichs, K., and Lewy, H., “On the partial difference equations of mathematical physics,” IBM J. Res. Dev. 11(2):215–234, 1967. 56. Yilmaz, E., Tian, T., Wong, V.W., and Heywood, J.B., “The contribution of different oil consumption sources to total oil consumption in a spark ignition engine,” SAE Trans. 1622– 1638, 2004. 197 57. Bilicki, Z. and Kestin, J., “Physical aspects of the relaxation model in two-phase flow,” Proc. R. Soc. London. A. Math. Phys. Sci. 428(1875):379–397, 1990. 58. Downar-Zapolski, P., Bilicki, Z., Bolle, L., and Franco, F., “The Non-Equilibrium Relaxation Model for One-Dimensional Flashing Liquid Flow, 3rd ASME,” JSME Joint Fluids Engineering Conference, 616, 1999. 59. Launder, B.E. and Spalding, D.B., “The numerical computation of turbulent flows,” Numerical prediction of flow, heat transfer, turbulence and combustion, Elsevier: 96–116, 1983. 60. Amsden, A.A. and Findley, M., “KIVA-3V: A block-structured KIVA program for engines with vertical or canted valves,” Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States), 1997. 61. Dursunkaya, Z., Keribar, R., and Richardson, D.E., “Experimental and numerical investigation of inter-ring gas pressures and blowby in a diesel engine,” SAE Technical Paper, ISBN 0148-7191, 1993. 62. Tian, T., Noordzij, L.B., Wong, V.W., and Heywood, J.B., “Modeling piston-ring dynamics, blowby, and ring-twist effects,” 1998. 63. Thirouard, B., “Characterization and modeling of the fundamental aspects of oil transport in the piston ring pack of internal combustion engines,” 2001. 64. Tamai, G., Experimental study of engine oil film thickness dependence on liner location, oil properties and operating conditions, 1995. 65. Noordzij, L.B., Measurement and analysis of piston inter-ring pressures and oil film thickness and their effects on engine oil consumption, 1996. 66. Petris, C. De, Giglio, V., and Police, G., “A mathematical model of the evaporation of the oil film deposed on the cylinder surface of IC engines,” SAE Technical Paper, ISBN 0148- 7191, 1997. 67. Tian, T., Rabute, R., Wong, V.W., and Heywood, J.B., “Effects of piston-ring dynamics on ring/groove wear and oil consumption in a diesel engine,” SAE Tech. Pap. 106(1997):1195– 1207, 1997, doi:10.4271/970835. 68. Inagaki, H., Saito, A., Murakami, M., and Konomi, T., “Development of two-dimensional oil film thickness distribution measuring system,” SAE Technical Paper, ISBN 0148-7191, 1995. 69. Nakashima, K., Ishihara, S., and Urano, K., “Influence of piston ring gaps on lubricating oil flow into the combustion chamber,” SAE Technical Paper, ISBN 0148-7191, 1995. 70. Saito, K., Igashira, T., and Nakada, M., “Analysis of oil consumption by observing oil behavior around piston ring using a glass cylinder engine,” SAE Trans. 1027–1038, 1989. 198 71. Furuhama, S., Hiruma, M., and Yoshida, H., “An Increase of Engine Oil Consumption at High Temperature of Piston and Cylinder,” SAE Trans. 3008–3017, 1981. 72. Orrin, D.S. and Coles, B.W., “Effects of engine oil composition on oil consumption,” SAE Trans. 475–490, 1971. 73. Stewart, R.M. and Selby, T.W., “The relationship between oil viscosity and engine performance-a literature search,” The Relationship Between Engine Oil Viscosity and Engine Performance, ASTM International, 1977. 74. Wahiduzzaman, S., Keribar, R., Dursunkaya, Z., and Kelley, F.A., “A Model for Evaporative Consumption of Lubricating Oil in Reciprocating Engines,” SAE Technical Paper, ISBN 0148-7191, 1992. 75. Audette III, W.E. and Wong, V.W., “A model for estimating oil vaporization from the cylinder liner as a contributing mechanism to engine oil consumption,” SAE Trans. 1080– 1092, 1999. 76. Zhang, Q., Tian, T., and Koeser, P., “A One-Line Correlation for Predicting Oil Vaporization from Liner for IC Engines,” SAE Tech. Pap. 2018-April:1–15, 2018, doi:10.4271/2018-01-0162. 77. Mittler, R., Mierbach, A., and Richardson, D., “Understanding the fundamentals of piston ring axial motion and twist and the effects on blow-by,” ASME 2009 Internal Combustion Engine Division Spring Technical Conference, American Society of Mechanical Engineers Digital Collection: 721–735, 2009. 78. Liu, L., Tian, T., and Rabuté, R., “Development and applications of an analytical tool for piston ring design,” SAE Technical Paper, ISBN 0148-7191, 2003. 79. Cheng, C., Kharazmi, A., and Schock, H., “Modeling of Piston Ring-Cylinder Bore-Piston Groove Contact,” SAE Technical Paper, ISBN 0148-7191, 2015. 80. Dowson, D., Economou, P.N., Ruddy, B.L., Strachan, P.J., and Baker, A.J.S., “Piston ring lubrication—Part II. Theoretical analysis of a single ring and a complete ring pack,” Energy Conserv. through Fluid Film Lubr. Technol. Front. Res. Des. 23–52, 1979. 81. Keribar, R., Dursunkaya, Z., and Flemming, M.F., “An integrated model of ring pack performance,” J. Eng. Gas Turbines Power 113(382), 1991. 82. Tian, T., Noordzij, Lb., Wong, V.W., and Heywood, J.B., “Modeling piston-ring dynamics, blowby, and ring-twist effects,” J. Eng. Gas Turbines Power Vol. 120, 1998. 83. Chowdhury, S.S., Kharazmi, A., Atis, C., and Schock, H., “Three-Dimensional Multi-phase Physics-Based Modeling Methodology to Study Engine Cylinder-kit Assembly Tribology and Design Considerations-Part I,” SAE Technical Paper, ISBN 0148-7191, 2020. 84. Cheng, C., Schock, H., and Richardson, D., “The dynamics of second ring flutter and 199 collapse in modern diesel engines,” J. Eng. Gas Turbines Power 137(11), 2015. 85. Furuhama, S., Hiruma, M., and Tsuzita, M., “Piston ring motion and its influence on engine tribology,” SAE Trans. 2929–2941, 1979. 86. Alice, S.E. and David, C.W., Penalty functions Handbook of Evolutionary Computation, Section C 5.2, 1996. 87. Ali, Z., Dhanasekaran, P.C., Tucker, P.G., Watson, R., and Shahpar, S., “Optimal multi- block mesh generation for CFD,” Int. J. Comut. Fluid Dyn. 31(4–5):195–213, 2017. 88. Allahyari, M., Esfahanian, V., and Yousefi, K., “The Effects of Grid Accuracy on Flow Simulations: A Numerical Assessment,” Fluids 5(3):110, 2020. 89. Lucchini, T., Torre, A. Della, D’Errico, G., Montenegro, G., Fiocco, M., and Maghbouli, A., “Automatic mesh generation for CFD simulations of direct-injection engines,” SAE Technical Paper, ISBN 0148-7191, 2015. 90. Heywood, J.B., “Internal combustion engine fundamentals,” Pub McGraw Hill Int. Ed., 1988. 91. Krause, W., Spies, K.H., Bell, L.E., and Ebert, F., “Oil Separation in Crankcase Ventilation- New Concepts Through System Analysis and Measurements,” SAE Trans. 1577–1588, 1995. 92. Shore, P.R., “Advances in the use of Tritium as a Radiotracer for Oil Consumption Measurement,” SAE Trans. 552–563, 1988. 93. Froelund, K., Menezes, L.A., Johnson, H.R., and Rein, W.O., “Real-time transient and steady-state measurement of oil consumption for several production SI-engines,” SAE Technical Paper, ISBN 0148-7191, 2001. 94. Ariga, S., Sui, P.C., Bailey, B., Kumakiri, T., Osumi, Y., and Sakamoto, A., “On-line oil consumption measurement and characterization of an automotive gasoline engine by SO2 method,” SAE Technical Paper, ISBN 0148-7191, 1992. 95. Apple, J., Gladis, D., Watts, W., and Kittelson, D., “Measuring diesel ash emissions and estimating lube oil consumption using a high temperature oxidation method,” SAE Int. J. Fuels Lubr. 2(1):850–859, 2009. 96. Delvigne, T., “Oil consumption sources in a modern gasoline engine including contribution of blow-by separator and turbocharger: an experimental study based on the use of radiotracers,” SAE Int. J. Fuels Lubr. 3(2):916–924, 2010. 97. Soejima, M., Smith, E.H., Sherrington, I., and Wakuri, Y., “A review of solutions for the mechanism of oil consumption in internal combustion engines,” SAE Technical Paper, ISBN 0148-7191, 2007. 200 98. Wong, V.W. and Hoult, D.P., “Experimental survey of lubricant-film characteristics and oil consumption in a small diesel engine,” SAE Trans. 282–295, 1991. 99. Ariga, S., Sui, P.C., and Shahed, S.M., “Instantaneous unburned oil consumption measurement in a diesel engine using SO2 tracer technique,” SAE Technical Paper, ISBN 0148-7191, 1992. 100. Franklin, A., “Calibration,” Can that be Right?, Springer: 237–272, 1999. 101. Chong, E.K.P. and Zak, S.H., “An introduction to optimization,” John Wiley & Sons, ISBN 0471654000, 2004. 102. HEEDS Version 19.1 Users’ Manual, 2020, No Title. 103. Chase, N., Redemacher, M., Goodman, E., Averill, R., and Sidhu, R., “A benchmark study of optimization search algorithms,” Red Cedar Technol. MI, USA 1–15, 2010. 104. Knudsen, D.C. and Fotheringham, A.S., “Matrix comparison, goodness-of-fit, and spatial interaction modeling,” Int. Reg. Sci. Rev. 10(2):127–147, 1986. 105. Chowdhury, S.S., Schock, H., and Kharazmi, A., “The Effect of Ring-Groove Geometry on Engine Cylinder-Kit Assembly Using Three-Dimensional Multiphase Physics-Based Modeling Methodology-Part II,” SAE Technical Paper, ISBN 0148-7191, 2021. 106. Cheng, C., “Numerical modeling of the power cylinder system for internal combustion engine with an emphasis on ring pack design,” Michigan State University, ISBN 1321425368, 2014. 107. Fischer, P., Min, M., Rathnayake, T., Dutta, S., Kolev, T., Dobrev, V., Camier, J.-S., Kronbichler, M., Warburton, T., and Świrydowicz, K., “Scalability of high-performance PDE solvers,” Int. J. High Perform. Comput. Appl. 34(5):562–586, 2020. 108. Tas, M.O., Banerji, A., Lou, M., Lukitsch, M.J., and Alpas, A.T., “Roles of mirror-like surface finish and DLC coated piston rings on increasing scuffing resistance of cast iron cylinder liners,” Wear 376–377:1558–1569, 2017, doi:10.1016/j.wear.2017.01.110. 109. Kumar, V., Sinha, S.K., and Agarwal, A.K., “Wear evaluation of engine piston rings coated with dual layer hard and soft coatings,” J. Tribol. 141(3), 2019. 110. Vinoth, I.S., Detwal, S., Umasankar, V., and Sarma, A., “Tribological studies of automotive piston ring by diamond-like carbon coating,” Tribol. - Mater. Surfaces Interfaces 13(1):31– 38, 2019, doi:10.1080/17515831.2019.1569852. 111. Friedrich, C., Berg, G., Broszeit, E., Rick, F., and Holland, J., “PVD CrxN coatings for tribological application on piston rings,” Surf. Coatings Technol. 97(1–3):661–668, 1997. 112. Pandey, S.M., Murtaza, Q., and Gupta, K., “Friction and sliding wear characterization of ion chrome coating,” SAE Tech. Pap. 1, 2014, doi:10.4271/2014-01-0946. 201 113. Kamo, L., Saad, P., Bryzik, W., and Mekari, M., “Ceramic coated piston rings for internal combustion engines,” Proc. World Tribol. Congr. III - 2005 607–608, 2005, doi:10.1115/wtc2005-64343. 114. Rozario, A., Baumann, C., and Shah, R., “The Influence of a Piston Ring Coating on the Wear and Friction Generated during Linear Oscillation,” Lubricants 7(1):8, 2019. 115. Hoppe, S. and Kantola, T., “DuroGlide® - New Generation Piston Ring Coating for Fuel- Efficient Commercial Vehicle Engines,” SAE Tech. Pap. 2014-Janua, 2014, doi:10.4271/2014-01-2323. 116. Wan, S., Wang, H., Xia, Y., Tieu, A.K., Tran, B.H., Zhu, H., Zhang, G., and Qiang zhu, “Investigating the corrosion-fatigue wear on CrN coated piston rings from laboratory wear tests and field trial studies,” Wear 432–433(June):202940, 2019, doi:10.1016/j.wear.2019.202940. 117. Johansson, S., Frennfelt, C., Killinger, A., Nilsson, P.H., Ohlsson, R., and Rosén, B.-G., “Frictional evaluation of thermally sprayed coatings applied on the cylinder liner of a heavy duty diesel engine: pilot tribometer analysis and full scale engine test,” Wear 273(1):82–92, 2011. 118. Lou, M. and Alpas, A.T., “Characterization of Lubricated Friction Behavior of Thermal Spray Steel Coatings in Comparison with Grey Cast Iron,” Lubricants 8(1):9, 2020. 119. Blümm, M., Baberg, A., Dörnenburg, F.T.H., and Leitzmann, D., “Innovative Skirt Coatings for Gasoline and Diesel Engine Pistons,” MTZ Worldw. 77(2):50–55, 2016. 120. Demas, N.G., Erck, R.A., and Fenske, G.R., “Tribological evaluation of piston skirt/cylinder liner contact interfaces under boundary lubrication conditions,” Lubr. Sci. 22(3):73–87, 2010. 121. Demas, N.G., Timofeeva, E. V, Routbort, J.L., and Fenske, G.R., “Tribological effects of BN and MoS 2 nanoparticles added to polyalphaolefin oil in piston skirt/cylinder liner tests,” Tribol. Lett. 47(1):91–102, 2012. 122. Wang, Y., Brogan, K., and Tung, S.C., “Wear and scuffing characteristics of composite polymer and nickel/ceramic composite coated piston skirts against aluminum and cast iron cylinder bores,” Wear 250(1–12):706–717, 2001. 123. Gavrilov, K. V, Morozov, A. V, Seleznev, M. V, Rozhdestvenskii, Y. V, Khozeniuk, N.A., Doikin, A.A., and Hudyakov, V.S., “Evaluation of Anti-Friction Properties of Solid Lubricant Coatings for a Piston Skirt of a High-Force Diesel,” J. Frict. Wear 41(5):480– 485, 2020. 124. Westerfield, Z., Totaro, P., Kim, D., and Tian, T., “An experimental study of piston skirt roughness and profiles on piston friction using the floating liner engine,” SAE Technical Paper, ISBN 0148-7191, 2016. 202 125. Cho, D.-H., Lee, S.-A., and Lee, Y.-Z., “The effects of surface roughness and coatings on the tribological behavior of the surfaces of a piston skirt,” Tribol. Trans. 53(1):137–144, 2009. 126. Zabala, B., Igartua, A., Fernández, X., Priestner, C., Ofner, H., Knaus, O., Abramczuk, M., Tribotte, P., Girot, F., and Roman, E., “Friction and wear of a piston ring/cylinder liner at the top dead centre: Experimental study and modelling,” Tribol. Int. 106:23–33, 2017. 127. Buyukkaya, E. and Cerit, M., “Thermal analysis of a ceramic coating diesel engine piston using 3-D finite element method,” Surf. Coatings Technol. 202(2):398–402, 2007, doi:10.1016/j.surfcoat.2007.06.006. 128. Hildyard, R., Bewsher, S.R., Walker, J., Umer, J., Saremi-Yarahmadi, S., Pacella, M., Mohammadpour, M., and Offner, G., “Tribological enhancement of piston skirt conjunction using graphene-based coatings,” Proc. Inst. Mech. Eng. Part D J. Automob. Eng. 235(5):1330–1350, 2021. 129. Fox, I.E., “Numerical evaluation of the potential for fuel economy improvement due to boundary friction reduction within heavy-duty diesel engines,” Tribol. Int. 38(3):265–275, 2005. 130. Perchanok, M., “Modeling of piston-cylinder lubrication with a flexible skirt and cylinder wall,” SAE Trans. 2339–2354, 2000. 131. Jian, Z., Zhong-yu, P., Shi-ying, L., Sheng-wei, S., and Li-jun, D., “Investigation of wear behavior of graphite coating on aluminum piston skirt of automobile engine,” Eng. Fail. Anal. 97(December 2018):408–415, 2019, doi:10.1016/j.engfailanal.2019.01.012. 132. Zhang, J., Piao, Z., and Liu, S., “Influence of skirt profile structure of gasoline engine piston on the friction and wear characteristics under standard conditions,” J. Tribol. 140(2), 2018. 133. Meng, Z., Zhang, L., and Tian, T., “Study of break-in process and its effects on piston skirt lubrication in internal combustion engines,” Lubricants 7(11):98, 2019. 134. Suman, A., “Abradable dry powder coatings on piston assembly components,” Google Patents, US 2007007 1990A1; United States Patent and Trademark Office; Mar. 29; 2007. 135. Suman, A. and Shamis, D.A., “Clearance control coatings - Low cost, abradable, lubricious,” SAE Tech. Pap., 2010, doi:10.4271/2010-32-0077. 136. Digital Metrology Solutions (2021). OmniSurf3D (Version 1.04), No Title. 137. Shi, R., Wang, B., Yan, Z., Wang, Z., and Dong, L., “Effect of surface topography parameters on friction and wear of random rough surface,” Materials (Basel). 12(17):2762, 2019. 138. Mathia, T.G., Pawlus, P., and Wieczorowski, M., “Recent trends in surface metrology,” Wear 271(3–4):494–508, 2011. 203 139. Sayles, R.S., “How two-and three-dimensional surface metrology data are being used to improve the tribological performance and life of some common machine elements,” Tribol. Int. 34(5):299–305, 2001. 140. Cohen, D 2008, Michigan Metrology Glossary of Texture Parameters, accessed December,2020. 141. Tomanik, E., “Modelling of the asperity contact area on actual 3D surfaces,” SAE Technical Paper, ISBN 0148-7191, 2005. 142. Temam, R., “Navier-Stokes equations: theory and numerical analysis,” American Mathematical Soc., ISBN 0821827375, 2001. 143. Pan, J., Nigro, R., and Matsuo, E., “3-D modeling of heat transfer in diesel engine piston cooling galleries,” SAE Technical Paper, ISBN 0148-7191, 2005. 144. Guyan, R.J., “Reduction of stiffness and mass matrices,” AIAA J. 3(2):380, 1965. 145. Panayi, A.P., “Numerical models for the assessment of the cylinder-kit performance of four- stroke internal combustion engines,” 2009. 146. Patir, N. and Cheng, H.S., “Application of average flow model to lubrication between rough sliding surfaces,” 1979. 147. Tripp, J.H., “Surface roughness effects in hydrodynamic lubrication: the flow factor method,” 1983. 148. Zhu, D., Cheng, H.S., Arai, T., and Hamai, K., “A numerical analysis for piston skirts in mixed lubrication—Part I: Basic modeling,” 1992. 149. Zhu, D., Hu, Y.-Z., Cheng, H.S., Arai, T., and Hamai, K., “A numerical analysis for piston skirts in mixed lubrication: part II—deformation considerations,” 1993. 204