Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Journal bearing optimization and analysis using streamline upwind Petrov-Galerkin finite element method Zengeya, Miles 2009

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2009_fall_zengeya_miles.pdf [ 3.29MB ]
Metadata
JSON: 24-1.0067254.json
JSON-LD: 24-1.0067254-ld.json
RDF/XML (Pretty): 24-1.0067254-rdf.xml
RDF/JSON: 24-1.0067254-rdf.json
Turtle: 24-1.0067254-turtle.txt
N-Triples: 24-1.0067254-rdf-ntriples.txt
Original Record: 24-1.0067254-source.json
Full Text
24-1.0067254-fulltext.txt
Citation
24-1.0067254.ris

Full Text

Journal Bearing Optimization And Analysis Using Streamline Upwind Petrov-Galerkin Finite Element Method  by  Miles Zengeya B.S.M.E University of Arkansas, U.S.A., 1988 M.Sc., Loughborough University of Tech., England, 1992  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENT FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Mechanical Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) May 2009 © Miles Nyamana Zengeya, 2009  Abstract A three-dimensional finite element thermo-hydrodynamic lubrication model that couples the Reynolds and energy equations is developed. The model uses the streamline upwind Petrov-Galerkin (SUPG) method. Model results indicate that the peak temperature location in slider bearing is on the mid-plane well as when pressure boundary conditions are altered in such a way that the inlet/outlet pressure is higher than the side pressure. The adiabatic temperature profiles of an infinite and square sliders are compared. The wider slider shows a higher peak temperature. Side flow plays a major role in determining the value and position of the peak temperature. Model results also indicate peak side flow at a width-to-length ratio of 2. A method of optimizing leakage, the Flow Gradient Method, is proposed. The SUPG finite element method shows rapid convergence for slider and plain journal bearings and requires no special treatment for backflow in slider bearings or special boundary conditions for heat transfer in the rupture zone of journal bearings. A template for modeling thermo-hydrodynamic lubrication in journal bearings is presented. The model is validated using experimental and analytical data in the literature. Maximum deviation from measured temperatures is shown to be within 40 per cent. The model needs no special treatment of boundary conditions in the rupture zone and shows rapid and robust convergence which makes it quite suitable for use in design optimization models and in obtaining closed relations for critical parameters in the design of journal and slider bearings. Empirically derived simulation models for temperature increase; leakage; and power loss are proposed and validated using the developed finite element model and experimental results from literature. Predictions of temperature increase, leakage, and power loss are  ii  better than those obtained for available relations in the literature. The derived simulation models include five important design variables namely the radial clearance, length to diameter ration, fluid viscosity, supply pressure and groove position. The derived model is used to minimize a multi-objective function using weight/scaling factors and Pareto optimal fronts. The latter method is recommended as preferable, and Pareto diagrams are presented for common bearing speeds. Including the groove location in the optimization model is shown to have a significant effect on the results. The lower bound of groove location appears to result in preferred power loss/side leakage values. Significant power loss savings may be realized with appropriate groove location.  iii  TABLE OF CONTENTS  Abstract…………...………………………………………………………………….……ii Table of Contents…...………………………………………………………….…………iv List of Tables………..………………………………………………………..………....viii List of Figures…………………………………………………………..………..……......x List of Symbols, Abbreviations and Parameters…...………………………….…….…...xv Acknowledgements…………………………………..……………………….…..……..xix  Chapter 1 Introduction……………………………..………………………..…….……1 1.1  1.2  Thermal Effects in Slider and Journal Bearings…..………………….….….…….1 1.1.1  Historical Background……………………..………………….…..………1  1.1.2  Thermal Effects and Boundary Conditions……..……….……….……….2  Literature Review………………………………………..………………………..3 1.2.1  Slider Bearings and Thermo-hydrodynamic Lubrication…...…..………...3  1.2.2 Journal Bearing THD Modeling………………………………..…………3 1.2.3  Convection-Dominated Flows and Upwinding………………..…….........5  1.2.4  Design Optimization…………..……………………………..…………....6  1.2.5  Journal Bearings and Optimization Models………………..…….….……7  1.3  Motivation For The Research……………………………………………………11  1.4  Structure of Research Work…………………………..…………..…….………..12  Chapter 2 Three-Dimensional SUPG FE Formulation………………...…….……....14 2.1  Introduction……………………………………….……………………………...14  2.2  Governing Equations……………………………………………...………….….14  iv  2.2.1  Pressure Computations……………………….………………………….14  2.2.2  Energy and Viscosity Variation……………….…………………………18  2.3  Finite Element Formulation with Streamline Upwind Petrov-Galerkin…...….....24  2.4  Modeling Procedure…………………………………………….……….………30  2.5  Chapter summary………………………………………………..………….……34  Chapter 3 Application of (SUPG) FE to Slider Bearings............................................35 3.1  Introduction………….…………………………………………...……….……...35  3.2  Lebeck Slider Bearing Data....……………………………………..….…....……35  3.3  3.2.1  Model Convergence Performance...…….…………………....…..….…...36  3.2.2  Predicted Pressure Values……………………...………..…..…..…....….38  3.2.3  Thermal Profile……………………...…………………..…….….…...…40  3.2.4  Leakage Flow Pattern – Flow Gradient Method………...…............…….46  Chapter Summary……………..………………………………..….…..…...……51  Chapter 4 Streamline Upwind Petrov-Galerkin FE and Journal Bearings…...…... 52 4.1  Introduction…………………………….………………………………………...52  4.2  Model Formulation……………………..…………………………..…….….…..52  4.3  Simulation Results………………………….…………………...………….…....55 4.3.1  Mesh Density, Convergence and Performance………..…….…………...55  4.3.2  Pressure and Pressure Gradient Predictions…………..…….…....………59  4.3.3  Pressure Gradient and Effect on Dissipation…………..…...……….…...60  4.3.4  Dissipation Source Terms…………………...………..…………….…....63  4.3.5  Thermal Profiles……………………………...…………...……….……..66  v  4.3.5.1 Standard Error and Goodness of Fit……………………………………...68  4.4  4.3.6  Pressure Simulations of Costa’s Data………………........………………77  4.3.7  Flow Predictions…………...…………..……………...……………...….79  Chapter Summary and Conclusions…………………….……..…………..….…79  Chapter 5 Formulation of Empirically Derived Model and Simulation Results......81 5.1  Introduction…………………………………...……………………………….....81  5.2  Governing Equations………………….…………………………………......…..86  5.3  5.2.1  Temperature Increase………………………………………………….....86  5.2.2  Side Leakage…………………………………………………………......88  5.2.3  Power Loss Determination…………………………………………...…..89  5.2.4  Pressure Computations…………………………………………….……..90  5.2.5  Formulation of the Optimization Problem…………………………...…..91  5.2.6  Implementation in Matlab…………………………………………...…...92  Validation of Proposed Equations……………………………………...….…….94 5.3.1 Effect of Supply Pressure ps………………………………...………...…94 5.3.2  Temperature Increase ΔT(X)……………………………………………..96  5.3.3 Leakage Flow QL Predictions…………………………….…………….102  5.4  5.3.4  Variation of C on Performance Parameters ΔT and QL….……………..104  5.3.5  Temperature Increase ΔT, Side Leakage QL as Function of Diameter....107  5.3.6  Simulation Results with ΔT (f1) and QL (f2) as Objective Functions…...108  Power Loss PL (F1) and Leakage QL (F2) as Objective Functions…………...…115 5.4.1  Power Loss Predictions from Model……………………………...….....115  5.4.2  Determination of Weighting Parameters and Scaling Factors…...…..…116 vi  5.4.3 5.5  5.6  5.7  Pareto Optimal Front Method………………………….…..……….….117  Extended Optimization Model……………………………………..…….…….125 5.5.1  Leakage Flow Equation and Validation……………………………...…125  5.5.2  Power Loss Equation and Validation……………………………….…..131  5.5.3  Optimization Simulations……………………………………..………..134  Interpretation of Optimization Results………………………………..………..140 and X T  5.6.1  Pareto Optimal Fronts with X  5.6.2  Uncertainty in the Results………………...…………………………….141  , , ,  , , ,  ,  140  Chapter Summary………………………………..…………………....….….…141  Chapter 6 Summary, Conclusions and Future Work…………….…..…….….…..143 6.1  Summary and Conclusions………………………….…………..….….….…...143  6.2  Recommendations for Future Work………………………….…….…..……..146  Bibliography……………………………………………………..……….…….……..144 Appendices……………………………………………………….…….…..….….…..157 Appendix A  Velocity fields u, v and w………………..…………………….157  Appendix B  Summary of Martin’s Procedure.………..……….……………159  Appendix C  Optimization Procedure………………..………………………160  Appendix D  Tables….……………………………….………………………166  Appendix E  Verification Examples……………….…………………………174  vii  List of Tables Table 3.1  Input values for the model……………………………...……….……….36  Table 3.2  Mesh density and effect on prediction………………..…….…………....36  Table 3.3  Convergence criteria…..……………………………..………….…….…37  Table 3.4  Isothermal peak pressure and load-carrying capacity, k=0.4..…….........38  Table 3.5  THD peak pressure and load capacity comparisons………….…….…....40  Table 3.6  Dimensionless flow rate as a function of L/B……………..…..….…......47  Table 4.1  Convergence of Peak Pressure Predictions…………….………..……….56  Table 4.2  Computation Times for one Iteration…………….……………..………..57  Table 4.3  Model predictions of analytical and Dowson’s data….……….…….…...59  Table 4.4  Effect of pressure gradients ∂p/∂x, and ∂p/∂z on temperature rise ΔT…...59  Table 4.5  Standard Error and r2 correlation for Temperature Predictions……..…...69  Table 4.6  Predicted versus measured temperatures, Costa’s data……….…..……..75  Table 4.4  Predicted versus measured temperatures as function of supply press...…76  Table 5.1  Computation Times for Optimization Model Time Stepping…….…..….84  Table 5.2  Values a, b, and c computation of maximum temperature Tmax….....…...87  Table 5.3  Standard error and square of the Pearson correlation……………..……..95  Table 5.4  Local Reynolds number from a sample of the data…………...….….…..97  Table 5.5  Summary of ΔT(X) and QL data comparisons………………..…….…...104  Table 5.6 Table 5.7 Table 5.8  Typical Pareto Table for Bearing at 4000 rpm, W=10kN…….…...…..125 Comparison of variance between model and Dowson’s predictions…...130 Variation of Performance Parameters PL and QL………………….…..140  viii  Table D.1  Input data from Dowson and Song bearings………………….….......…166  Table D.2  Input values for the bearing tested by Boncompain et al. ……….....….169  Table D.3  Input values for Costa’s bearing………………………………...……...170  Table D.4  Input values for Ferron’s bearing………………………….……...……171  Table D.5  Additional input data to run optimization simulations…….…..…….…172  Table D.6  Input data for Dowson and Claro bearings………………….…..……...173  Table E1  Minimum and minimum values of the solution………………...………177  Table E2  Comparison of SUPG and Galerkin predictions………….…...………..179  Table E3  Prediction comparison, analytical vs SUPG FE………………...…...…182  Table E4  Input parameters and predictions……………………………...….…….183  Table E5  Input values for Mitsui bearing……………………………...……...…..184  Table E6  Maximum and average error for Mitsui’s data……………..………….185  Table E7  Input parameters for journal tested.…………………………..….…….186  Table E8  Error on thermal predictions, SUG vs THD model…………...………..189  Table E9  Performance parater predictions for Pierre’s data……..…………....…189  ix  List of Figures  Figure 1.1  Flow chart of a simple genetic algorithm……………………..……....…10  Figure 1.2  Thesis roadmap…………………………………………….…….………13  Figure 2.1a  Slider bearing………………………….……………………….………...16  Figure 2.1b  Journal bearing and notation………………………………….….…....…16  Figure 2.2  Inlet temperature variation for Dowson’s bearing data……….…......…..21  Figure 2.3  Three dimensional slider bearing geometry……………………...……....25  Figure 2.4a  Solution domain Ω showing boundary Γ1 and unit normal vector n…….26  Figure 2.4b  Transformation of domain Ω……………………………………....….....26  Figure 2.5  Flowchart for thermo-hydrodynamic lubrication model…………......…..33  Figure 3.1  Convergence of load capacity solution…………………………………..37  Figure 3.2  Non-dimensional isothermal pressure predictions……….……....………39  Figure 3.3a  Normalized temperatures for 2-D SUPG model………….…..….………42  Figure 3.3b  Temperature predictions for 3-D SUPG model………………..…......….42  Figure 3.4  Contour map of thermal field, L/B=10….…………………...…………..43  Figure 3.5  Thermal field with L/B=0.5……...………………...………..…..….…....44  Figure 3.6a  Model isotherms for Lebeck’s data……………………………...….........45  Figure 3.6b  Comparison of temperature profiles from SUPG and Kumar………..….45  Figure 3.7  Dimensionless side leakage QL as a function of L/B, k=0.4………......…48  Figure 3.8  Flow Gradient, k=0.4…………………………………...………….….....48  Figure 3.9  Side leakage as a function of L/B, k=0.4……………….……..……...….49  Figure 3.10  Load capacity as a function of k……………………………..…….….....50 x  Figure 4.1a  Schematic representation of journal bearing……..……………………....53  Figure 4.1b  Journal bearing unwrapped…………………………...………………….53  Figure 4.1c  Isometric view of finite element mesh………………...…………………54  Figure 4.2  Cross section of the film region……………………………………...…..54  Figure 4.3  Convergence analysis on Dowson’s data…………………….…..…...….56  Figure 4.4  Peak Pressure Prediction for Dowson’s exp. Value of 35MPa…………..57  Figure 4.5  Model, experimental, and analytical predictions, Dowson’s data……….58  Figure 4.6  Convergence analysis on pressure gradient dp/dz……………………….61  Figure 4.7a  Dimensionless pressure predictions at bearing mid-plane, model.............62  Figure 4.7b  Dimensionless pressure predictions at bearing mid-plane, Bonc............. 62  Figure 4.7c  Pressure map showing the rupture and reformation interfaces………..…62  Figure 4.8a  Pressure gradient ∂p/∂x variation along the bearing z axis…………..…..63  Figure 4.8b  Pressure gradients ∂p/∂x, ∂p/∂z variation along the bearing x axis……...63  Figure 4.9a  Velocity gradient ∂u/∂y as a function of X …………………………….65  Figure 4.9b  Magnified view at the peak gradients……………..…….………….…....65  Figure 4.10  Velocity gradient ∂w/∂y as a function of X ……………………..…..….66  Figure 4.11a Thermal profile, Boncompain’s data………………..…….……………..71 Figure 4.11b Figure 7 from Boncompain et al. ………………………………..…...…71 Figure 4.12a,b Thermal profile in the inlet, cavity and bearing for Costa’s data….…...72 Figure 4.13  Bush thermal profile……………………………………………...……...73  Figure 4.14  Thermal predictions at the bush/oil interface – Costa’s data………........74  Figure 4.15  Predicted versus measured peak temperatures………………..…………76  xi  Figure 4.16  Pressure prediction as a function of groove location θ………….…........78  Figure 4.17  Flow predictions: model, measured data and Martin’s formula…............79  Figure 5.1  Flowchart for optimization model……………………………..…..…….85  Figure 5.2  Model versus Khonsari’s predictions – Costa’s data…………..….....….95  Figure 5.3  ΔT(X) as function of speed Ns for Dowson’s bearing…………...………98  Figure 5.4  ΔT(X) as a function of speed Ns for Ferron’s bearing…………..….........99  Figure 5.5  ΔT(X) as function of speed Ns, Costa’s bearing…………………..….…100  Figure 5.6a  ΔT(X) as function of clearance ratio C, Dowson data…………..…...…101  Figure 5.6b  ΔT(X) as function of diameter D, Dowson data………………..…..…..101  Figure 5.7  Side leakage predictions for Dowson’s bearing………………....….......102  Figure 5.8  Side leakage predictions – Ferron’s bearing…………………..………..103  Figure 5.9  Side leakage predictions – Costa’s bearing…………………..….……..103  Figure 5.10  Side leakage as a function of clearance C, speed 25 to 201 rps…....…..105  Figure 5.11  ΔT(X) versus clearance ratio C, Dowson’s bearing, Ns=201 rps…….…106  Figure 5.12  Side leakage as a function of clearance and three different speeds….....106  Figure 5.13  ΔT(X) as a function of diameter D……………………………….….....107  Figure 5.14  QL as a function of diameter D, Dowson’s bearing, Ns=25 rps…...…...108  Figure 5.15a Objective function from Song et al. ……………………………..…....111 Figure 5.15b f(X) predictions from SQP and hybrid algorithm………………..….…111 Figure 5.16a Clearance ratio C as a function of speed for Song’s model……..….....112 Figure 5.16b Clearance ratio C as a function of speed, current model…………....…112 Figure 5.17a Length to diameter ratio λ as function of speed for Song’s model…....113  xii  Figure 5.17b Length to diameter ratio λ as function of speed, current models…….....113 Figure 5.18a Viscosity variation, Song’s model……………………………………...114 Figure 5.18b Viscosity variation versus speed, hybrid and SQP………….…..…...…114 Figure 5.19  Power loss prediction – Dowson’s data……………………...…......…..115  Figure 5.20  Power loss predictions – Costa’s data………………………...…...……116  Figure 5.21  f(X) as function of α1/α2, scaling parameters: β1=1, β2=105 …..…........118  Figure 5.22  Model predictions with α1/α2=1/5, β1=1/400, β2=105 ………....……...119  Figure 5.23  Model predictions with α1/α2=1/5, β1=1/35, β2=105 ……………….....120  Figure 5.24  Model predictions with α1/α2=1/5, β1=1, β2=105 ……………...…...…121  Figure 5.25  Pareto optimal front for speeds of (a)1500 rpm and (b) 2000 rpm…......123  Figure 5.26  Pareto optimal front for (a) 3000 rpm and (b) 4000 rpm…...…….....….124  Figure 5.27a Leakage flow prediction and data as function of supply pressure……...128 Figure 5.27b Data from Costa et al. ……………………………………………….....128 Figure 5.28  Leakage flow as function of supply pressure, Costa………...……….…129  Figure 5.29  Non-dimensional leakage flow as function of eccentricity ratio…….…129  Figure 5.30  Leakage flow as function of speed, ps0= and 170 kPa, Dowson…….…130  Figure 5.31  Non-dimensional leakage flow as function of eccentricity ratio….……131  Figure 5.32  Torque predictions as function of groove position, Costa……..….……133  Figure 5.33  Torque predictions as function of groove non-dim pressure……….…..133  Figure 5.34  Variation of torque as function of speed, Costa…………...…..…..……134  Figure 5.35  Variation of XT as a function of speed, W=10kN…..……………..……135  Figure 5.36  Variation of XT as a function of speed, W=700N…..……………..……138  Figure 5.37a,b Pareto optimal front for speed of 1500 rpm, W=10kN….………..…….139  xiii  Figure D1  Thermocouple location for pressure and temp. measurements…….…...168  Figure E1  Region of the skew…………………………………………………..….176  Figure E2  Mesh for convection-diffusion problem………………………….…….177  Figure E3  Three-dimensional plots…………………………………………..…….178  Figure E4  SUPG and Galerkin predictions as a function of speed………….……..180  Figure E5  Slider geometry ……………………………………………...….……...181  Figure E6  SUPG predictions versus measured temperatures, Mitsui………….…..185  Figure E7  Experimental and THD temperature predictions as function of angle....187  Figure E8  Comparison of SUPG and Pierre’s THD model………………………..188  xiv  List of Symbols, Abbreviations and Parameters B length of slider (x-dimension) Bib, Bis Biot number for bush, shaft cf, cb, cs specific heats of lubricant, bush and shaft, [J/kg.K] C radial clearance, [m] D journal diameter, [m] πD Bearing circumferential length [m] e  eccentricity, [m]  Ec Eckert number Ec = U2/(cfΔT) F objective function;  Ffr frictional force, [N]  f coefficient of friction;  fR reaction load on bearing, [N]  h film thickness, [m];  hp thickness of the pad  Hb,s convective heat transfer coefficients for bush, shaft, [W/m2.k] IP point of inflection where gradient is zero k ratio of slider inlet (hi) to outlet (ho) dimensions (k=hi/ho) kf, kb, ks, kg thermal conductivities fluid, bush (pad), shaft (runner) and gas, [W/m.K] L width of bearing perpendicular to the direction of motion (z-dimension), [m] n outward normal vector N mesh density (number of elements in the mesh) Ni the ith basis (shape) function Nu Nusselt number, Nu=(HbR/kb) Nv number of volume (brick) elements Ns journal speed, rev/s  xv  Pe Pecklet number, Pe = ρocpωC2/kf PL Power loss, [W] p film pressure, [N/m2];  po ambient pressure [Pa];  ps supply pressure [Pa]  Pr Prandtl number, Pr = cpμ/kf Qs leakage flow rate [m3/s];  Qv the volume flow rate [m3/s]  R1, R2, R3 bush inner, shaft, and bush outer radii respectively, [m] (see figure 2.1b) r2 the square of the Pearson product moment correlation coefficient T lubricant temperature, [K];  Ta, Ti ambient, inlet temperature, [K]  Ts, Tb, Tr shaft, bush and runner temperature, [K] Tfr frictional torque [N.m] U shaft surface velocity, [m/s] u velocity vector u, v, w fluid velocity component in x, y, z-axis respectively, [m/s] X vector of design variables W load capacity of the bearing, [N] x, y, z coordinates in the circumferential, radial, and axial bearing directions z1 parameter for computing attitude angle ϕ  α thermal diffusivity of the lubricant (α=kf/(ρcf)) α1, α2 Weigthing parameters, Optimization Model β temperature-viscosity coefficient [1/K] β1,β2 Scaling parameters, Optimization Model ε eccentricity ratio, ε = e/C εr target error magnitude xvi  ζ streamline upwind parameter μ fluid viscosity [Pa.s] δ lubricant thermal expansivity [1/K], κ1, κ2 temperature rise parameters λ, λs length to diameter ratio for journal (L/D) and slider (L/B) bearings respectively ρ density of lubricant [kg/m3] σd standard deviation of predictions θ circumferential coordinate of angle from oil groove position (Fig. 1), [radians] ϕ attitude angle (Fig. 1), [radians] ϕcav angle from oil groove to the cavitation boundary/interface, [radians] Ω  volume domain, discretized for finite element; Γ boundary to domain Ω  Subscript a  indicates ambient conditions  i  indicates conditions at the oil inlet  o  indicates conditions at the outlet of slider (trailing edge)  Superscript ‘¯’ indicates a dimensionless parameter (note: ‘+’ is used for normalized temperature  T+ =  T (in °K) which is a ratio to show the magnitude of temperature increase. Ti  Also μ+ as defined in non-dimensional parameters below) e denotes elemental values s denotes a shape function for surface element (finite element formulation)  xvii  se denotes number of integration points, surface element ve denotes number of integration points, volume element  Non-dimensional parameters (those in brackets are specific to sliders). Note that symbols are defined above u =  u ; U  v=  T+ =  T ; Ti  ρ+ =  ρ ; ρi  Pr E c =  Bib =  (Q =  μ iU 2 k f Ti  R v ; CU  Ts =  w=  w x ; x= ; U πD  H b R3 ; kb Q ); ⎛ ho ⎞ 2 B U⎜ ⎟ ⎝B⎠  y z h ; z= ; h = h L C  Ts ; T = β (T − Ti ) ; β = β Ti ; δ = δTi Ti  ρ = 1 − δ (T + − 1) ;  ;  y=  C2 p ; p= μR 2 N s  Bis =  μ = exp[− T ];  (p=  μ=  μ μ eff ;  μ+ =  pho2 ) (for sliders) μ iUB  HsL ks  ( Tr =  Tr ho ); μaUB2  (W =  W ⎛ B3 μ iU ⎜⎜ 2 ⎝ ho  ⎞ ⎟⎟ ⎠  )  For the bush and pad (see fig. 2.1b): ,  ,  ,  ,  For the shaft: where ys is measured from the origin as indicated in figure 2.1b. Other non-dimensional values presented equation 2.4 and 2.7.  xviii  μ μi  Acknowledgements  This thesis would not have been possible without the assistance, support and encouragement from those who interact with me on a regular basis.  Top on that list is my supervisor Professor Mohamed S. Gadala for his academic and financial support. His guidance and encouragement have been crucial in finishing this work.  Special mention goes to Dr Guus Segal whose sharp numerical modeling abilities were invaluable in implementing the model in the finite element method. His proof reading of the journal papers is much appreciated. Sabrina Crispo, a co-worker at Koerner’s Pub who helped in proof reading the papers.  I also would like to thank Professor Michel Fillon, Université de Poitiers, France for kindly agreeing to include figures from some of the papers he contributed.  Finally I am grateful to my family whose moral support at crucial times made a difference in this work. And yes dad, you made such a difference in my life. You continue to inspire me even when times are hard and quitting seems the only option.  Miles Zengeya  xix  Chapter 1 Introduction and Literature Review  1.1  Thermal Effects in Slider and Journal Bearings  Slider and journal bearings are widely used in applications such as mechanical seals, collar thrust bearings, machine tool ways, automotive engine crankshaft and connectingrods, steam turbines in power generating stations, railway wheels, piston rings, and many others. They have high reliability, good load-carrying capacity, excellent stability and durability and low coefficient of friction. Relative motion between two surfaces leads to not only an increase in temperature due to friction and viscous energy dissipation but thermal gradients that can develop on the surfaces. Almost all the frictional energy is dissipated in the form of heat. A thorough understanding of thermal effects in bearings is important in order to predict bearing performance. Investigation of the lubricant temperature fields leads to a better understanding of the load-carrying capacity of bearings and how complete loss of bearing clearance known as seizure can be avoided [1].  1.1.1  Historical Background  Historically hydrodynamic lubrication theory was developed without consideration of thermal effects [2]. Isothermal models were used and adjusted for thermal effects and these models were adequate for relatively slow machinery where viscous energy was not very high. However, with modern machinery that runs at high speeds and loads thermal effects become significant. Not only must the non-linear effects associated with material property changes be addressed, the system thermal effects which couple the fluid film  1  with the overall mechanical system and its environment must be taken into account. The thermo-hydrodynamic (THD) problem is not easy to deal with considering the difficulty in extracting from the plethora of variables and boundary conditions those few parameters which embody the essentials of the process. Another problem is the difficulty in manipulating the complex set of governing differential equations – the Reynolds and energy equations, and other auxiliary material property expressions as an interconnected and coupled system. Lastly, the problem must be formulated such that solutions obtained are applicable to as many bearing configurations as possible.  1.1.2 Thermal Effects and Boundary Conditions Temperature variations are due to several causes, one of which is the heat generated by viscous dissipation. Temperature increase due to this phenomenon leads to a change in fluid properties, particularly its viscosity which lowers which in turn leads to lower loadcarrying capacity of the film; this might cause the two surfaces to come into contact and eventually cause bearing seizure. A properly designed bearing will operate in the range where fluid property variations are in a range that will not cause such severe effects. The second mechanism is heat transfer to and from the fluid film which depends on material properties of the shaft and bush as well as the fixtures that keep the two in place. The shaft might be operating in high temperature environments that may result in heat transfer to the film, or the bush may be surrounded by hot ambient conditions that influence heat transfer at the interface. A third mechanism that affects the thermal map is the inlet mixing conditions. The hot recycled oil mixes with a fresh stream at the supply groove. Many mixing models appear in the literature with varying degrees of prediction accuracy [3,4,5,6].  2  Thermal aspects in fluid film bearings are important due to [2]: a) Their effect on bearing performance as mentioned above. This includes film thickness h, power loss PL, and side leakage rates QL. b) Magnitude and location of the maximum temperature. c) Effects of thermal gradients on the geometry of the fluid film (thermal distortion). d) Heat flow to and from the components of the system. Although their importance may not be the same, consideration of the above aspects plays some part in most thermo-hydrodynamic (THD) analysis. Besides determining performance and possible seizure due to degradation of the film, thermal gradients can also cause failure of the bearing material by instigating cracking. It is generally believed that proper design of journal and slider bearings is by appropriate consideration of thermal effects. Thermal effects, however, are not the only focus of the work carried out here. The bearing is viewed as a system and optimization of variables that designers have control over will also be investigated. Literature review of the state of thermohydrodynamic lubrication follows.  1.2  Literature Review  1.2.1  Slider Bearings and Thermo-hydrodynamic Lubrication  Many studies on slider bearings appear in the literature [1,2,4,7,8,9,10,13,14]. Recent examples include the work of Rodkiewicz [7], Schumack [8] and Kumar et al. [9,10]. Schumack derived a spectral element scheme for thermo-hydrodynamic (THD) lubrication whereas Kumar et al. [9] used two dimensional finite elements to numerically analyze the load-carrying capacity of slider bearings under thermal influence by solving the coupled nonlinear partial differential equations governing fluid film temperature field  3  with either isothermal or adiabatic boundary conditions. At high speeds convection dominates the flow resulting in numerical oscillations due to the convective terms in the energy equation. Kumar et al. [9] used a finite element numerical scheme with streamline upwind Petrov-Galerkin (SUPG) as proposed by Brooks and Hughes [11] and Grygiel [12] to formulate a two dimensional version. Zengeya and Gadala [13] used a control volume finite difference (CVFD) scheme that recognizes location of the nodes with respect to fluid flow and applies appropriate weighting to solve the 2D thermohydrodynamic problem. The models mentioned above mainly focus on the two dimensional solution of the coupled equations. This works well for predicting the behavior of infinitely wide bearings. One important aspect of the work to be considered is modeling bearings in three dimensions and investigating the effect of aspect ratios L/B and L/D on leakage flow. Both finite and infinite bearings would be analyzed. Although other workers used threedimensional models for slider bearings [14,15,16,17,18], the emphasis here would be on formulating the thermo-hydrodynamic problem in a way that applies to both sliders and plain journal bearings with little modification and to use the new template in deriving optimization equations for journal bearings.  1.2.2 Journal Bearing THD Modeling The field of application of journal bearings is large. Crankshafts and connecting-rods in automotive engines are two common examples; they are also used extensively in steam turbines and railway wheels. Properly installed and maintained, journal bearings have essentially infinite life. However, a thorough investigation of the lubricant temperature fields leads to a better understanding of the load-carrying capacity of the bearing and to  4  minimization of seizure. Higher temperatures affect bearing performance and reduce the margin against seizure. There are numerous studies on thermal effects in journal bearings. Khonsari [1] gave a good review of much of the work up to 1992. The work edited by Pinkus [2] is also valuable as a reference on this important area. Numerical solution of the thermal problem in journal bearings can be broadly classified into the finite volume (FV), finite difference (FD) and finite element (FE) methods. Significant contributions in the work done using finite volume is presented by Paranjpe and Goenka [19], Paranjpe and Han [20,21]. Stachowiak and Batchelor [22] present a numerical CVFD model while Booker and Huebner [23], Goenka [24], Labouff and Booker [25], Bayada et al. [26], and Kucinschi et al. [27] use finite elements. An important consideration in journal bearing modeling is cavitation, and some of the most significant contributions in this area were made by Elrod [28], Kumar and Booker [29,30], and more recently by Optasanu and Bonneau [31].  1.2.3  Convection-Dominated Flows and Upwinding  The lubricant flow problem is convection dominated, resulting in numerical oscillations or ‘wiggles’ if the symmetrically weighted Galerkin method is used. Some form of upwinding whereby nodes upstream of the flow direction are assigned a higher weighting than downstream ones is required. Gero and Ettles [16] solve the problem using a backward difference Galerkin scheme. Colynuck and Medley [32] use a control volume finite difference (FD) scheme that uses compass directions (east, west, north, south) to enforce flow directionality into the controlling FD equations. Kumar et al. [9] use a finite element numerical scheme with streamline upwind Petrov-Galerkin (SUPG) scheme for two dimensional THD problems.  5  It should be emphasized that there should be a way to define upstream and downstream positions in such algorithms. A method to account for possible backflow at the inlet must also be included.  1.2.4  Design Optimization  Tied to the thermal performance of journal bearings is the choice of variables a designer may control. A designer may choose the type of oil and its initial viscosity for both slider and journal bearings. Other parameters that are determined at the early stages of the design process include the choice of materials for the two surfaces, general bearing dimensions, the speed at which the bearing is expected to run, load and oil supply pressure.  Choosing these variables influences a second group of variables, usually  referred to as performance parameters group, that the designer has no direct control over, and these may include: the coefficient of friction f, temperature rise ΔT, volume flow rate Q (which in turn affects leakage flow (QL), and the minimum film thickness hmin. Certain limitations on the performance parameters must be imposed by the designer to ensure satisfactory performance. One can surmise the design of journal and slider bearings as choosing the first group of variables such that limits on the second group are not exceeded [33]. The second part of the thesis develops an optimization scheme to help designers select the best possible set of design variables for a particular application. In order to properly carry out the optimization procedure, certain relations of key parameters must be established. In the work of the second part of the thesis such relations for leakage rate and power loss are proposed, tested and verified. A brief review of literature on journal bearing optimization follows.  6  1.2.5  Journal Bearings and Optimization  Multi-objective optimization of bearings can be viewed as a minimization of temperature increase and side leakage as proposed by Hashimoto [34], Hashimoto and Matsumoto [35] and used extensively by Song et al. [36], and Yang et al. [37] for high-speed, short journal bearings (0.2<λ<0.6). Hashimoto uses weighting and scaling factors to combine the two objective functions into one multi-objective function. However this has the disadvantage that the designer must have prior knowledge about the relative importance of each objective. Hirani and Suh [38] propose minimizing the leakage and power loss since these two are ‘independent axioms’ [39]. Hirani [40,41] utilizes an alternative approach to multi-objective optimization, namely Pareto optimal fronts. This approach uses a posterior articulation of the weights in that the designer initially generates a number of non-inferior (a set of equally efficient) solutions from which a final decision is made on any one solution. These two methods can be formally defined as [42] •  Prior articulation of preferences: involves combining individual objective functions using pre-decided weight factors into a single utility function [34,35,38,43] and solving problem as a single objective optimization problem. Weight factors express relative importance of objective function in the overall utility measure.  •  Posteriori articulation of preferences: involves the optimizer generating a number of non-inferior (a set of equally efficient) solutions from which a final decision is made on any one solution. This approach is often referred to as Pareto optimal approach. As a set of many trade-off solutions are already available with their pros and cons, Pareto optimal approach helps high level qualitative decisions.  7  Many optimization solution techniques are available depending on the nature of the objective function [42]. A popular classical solution technique for multi-objective functions with nonlinear constraints is the sequential quadratic programming (SQP). The method uses a quasi-Newton updating procedure and approximates the Hessian matrix of the optimization problem by an identity matrix. Various techniques may be used to improve the calculation of the search direction and to control the step size. Also, various modifications may be used to improve the identity matrix assumption of the Hessian matrix. Based on the work of Powell [44] the method mimics Newton’s method for constrained optimization whereby at each major iteration an approximation is made of the Hessian which is based on knowledge from previous steps. The method requires an initial guess for the design variables. Genetic and evolutionary algorithms are the current state of the art solution techniques for optimization problems with discontinuous objective functions or those with many local minima which make SQP unsuitable [44,45]. The term evolutionary suggests the natural process associated with biological evolution, particularly the Darwinian rule of selection of the fittest individuals in a population. Genetic algorithms (GA) were first proposed by Holland [46] and extended further by DeJong [47] and Goldberg [45]. They are an efficient search technique which apply the rules of natural genetics to explore a given search space [48]. They are robust, adaptive and well behaving for problems with combination of complex, discontinuous and discrete functions. Genetic algorithms maintain a population of encoded solutions, and guide the population towards the optimum solution [47]. Thus they search the space of possible individuals and seek to find the best fitness strings. Rather than start from a single point within the search space  8  as in SQP, GA’s start with an initial set of random solutions of population. The solutions are represented by strings (chromosomes) which are coded as a series of zeros and ones or a vector solution. GA’s are non-deterministic and do not require the evaluation of objective function differentials, hence their classification as zero-order optimization methods [49] which require only function evaluation. Figure 1.1 shows a simple flowchart for a genetic algorithm [50]. The process starts with the selection of a randomly chosen population and evaluation of fitness values for the individuals. An iterative process follows until satisfaction of the termination criteria. Within the reproduction phase, each individual’s fitness is evaluated, genetic operators that include crossover and mutation are applied to produce the next generation. The new individuals replace the parent generation and re-evaluation of the fitness of new individuals performed. A new set of chromosomes are produced at each generation using information from the previous generation. Related to the two main solution techniques is a group of ‘hybrid schemes’ that combine the benefits of both SQP and GA to solve optimization problems. The genetic algorithm is used to identify the area of the global minimum and SQP finds the exact location of the minimum [51]. Hybrid schemes are popular for solving journal bearing problems where the objective function has local peaks and valleys or is discontinuous and use of SQP alone may lead to erroneous solutions.  9  SEEDING  Start  Input: Design variable coding Initial population Objective Function  EVALUATION  REPRODUCTION  Generation 1  Perform selection Parent I Parent II  Perform crossover Perform mutation  Until temporary population is full  Perform other genetic operators  Updating existing generation  Evaluation  Termination criteria satisfied  Yes  Write final result  No Generation = Generation + 1  End  Figure 1.1 Flow chart of a simple genetic algorithm [50].  10  1.3  Motivation For The Research  The THD models described in section 1.2.1-1.2.3 are generally specifically formulated for slider or plain journal bearings. This approach obviously works well for dedicated computations in slider or journal bearings. The need for a general thermo-hydrodynamic formulation applicable to both slider and plain journal bearings with minimum modification gives motivation for the work. A holistic approach is adopted whereby the theoretical equations are formulated, validated and used in analyzing bearing phenomenon. Encompassed in this approach is optimization of design variables that engineers can control. A formal statement of the research objectives: •  to formulate the three-dimensional finite element THD solution in a form applicable to both slider and journal bearings, is robust, and converges rapidly to be used in validating optimization equations developed in subsequent sections (chapter 5). Key elements of the model include wide applicability, to be used for slider as well as journal bearings, ease of applying boundary conditions, and ease of implementation in batch runs for optimization procedures.  •  The model should have minimal requirements for boundary conditions that are easy to implement. No need for complex boundary conditions to account for dissipation and thermal fields in the rupture zone as explained in section 4.3 where pressure gradients are analyzed in detail.  •  Is based on the streamline upwind Petrov-Galerkin finite element formulation for its accuracy in convection-dominated flows and its ability to account for backflow in the groove region.  11  •  The idea of a ‘thermo-hydrodynamic bearing template’ is proposed as a direct outcome of the general formulation proposed. The template is easy to implement and gives reasonable predictions for bearing operating characteristics as shown in subsequent chapters.  •  The template is implemented such that groove location θ can be changed which requires generating a new mesh. This novel dynamic aspect of generating output proved critical in validating the optimization model in the second part of the work.  •  Empirically derived bearing optimization equations are proposed as a major contribution of this work. The equations are derived from first principles and validated with the aid of the developed template.  1.4  Structure of Research Work  In the first part of the thesis, a robust thermo-hydrodynamic (THD) solution scheme applicable to both slider and plain journal bearing is proposed. The scheme, developed in the finite element method, has to be robust and exhibit rapid convergence characteristics due to the iterative nature of the optimization procedures proposed in the second part. Application of the scheme to slider bearing follows with analysis and discussion on certain aspects of obtained results. Chapter 4 shows how the scheme can be applied to plain journal bearings with minor adjustments. Application and discussion of results obtained for journal bearings are presented. The second major section of the work focuses on using the developed finite element scheme to formulate optimization model for slider and journal bearings. New relations for leakage rate and power loss are proposed, tested and verified. Data on journal  12  bearings from the literature is used to simulate and validate the optimization equations. Two main optimization solution methods are compared and recommendations on the preferred method made in the final section. Figure 1.2 summarizes the flow of the work carried out.  Develop robust THD scheme suitable for slider and plain journal bearings  Validate model with slider bearing data from literature  Extend THD model to plain journal bearings and validate with data from literature  Use developed THD scheme to develop and validate an Optimization procedure for plain journal bearings  Use developed Optimization Scheme to make Design Recommendations via Pareto Charts  Report on optimization results & conclusions  Figure 1.2 Thesis roadmap  13  Chapter 2 Three-Dimensional SUPG FE Formulation 2.1  Introduction  The previous chapter explored literature in thermo-hydrodynamic lubrication and optimization. This chapter develops the three-dimensional model that extends current models. A general scheme that can be used for both slider and journal bearings is developed. Governing equations for hydrodynamic lubrication and energy dissipation are initially presented before formulation of the streamline upwind Petrov-Galerkin finite element method. The formulation represents the first contribution in this work [52,53].  2.2  Governing Equations  2.2.1  Pressure Computations  The Reynolds equation describes pressure variations in thin film flow where the film thickness h is very small compared to the size of the bearing. One can show by dimensional analysis that the continuity and Navier-Stokes equations describing fluid flow in this case may be approximated by the Reynolds equation for incompressible flow. Assumptions made in deriving Reynolds equation are: 1. Pressure is invariant across the fluid film (in the y direction, see figure 2.1) 2. Inertial forces are negligible 3. The flow is laminar 4. The lubricant behaves as a Newtonian fluid and there is no slip at the boundaries. The Reynolds equation takes the form  14  ⎛ h3 ⎞ (∇p − γfb ) + h (U1 − U 2 ) ⎟⎟ + dh + k ( p − po ) = 0 ∇⎜⎜ − 2 ⎝ 12μ ⎠ dt  (2.1)  where ∇ denotes the divergence, h and dh/dt are the film thickness and rate of change of the film thickness, μ the fluid viscosity (Pa.s), u the velocity vector of the runner and pad (figure 2.1a), p and po the gauge and ambient pressure respectively, fb body forces, and γ and k are constants. Figure 2.1 illustrates a slider and journal bearing with notations. Dowson [54] derived a generalized form of the Reynolds equation which allows variations in fluid properties along and across the film. Dowson further simplified the momentum equation by an order of magnitude analysis and neglecting body forces  F1 ⎞ F1 ⎤ ∂ ⎧ ∂p ⎫ ∂ ⎧ ∂p ⎫ ∂ ⎡ ⎛ ∂h + (v 2 − v1 ) (2.2) ⎨F2 ⎬ + ⎨ F2 ⎬ = ⎢U 2 ⎜⎜ h − ⎟⎟ + U 1 ⎥ − U 2 ∂x ⎩ ∂x ⎭ ∂z ⎩ ∂z ⎭ ∂x ⎣ ⎝ F0 ⎦ F0 ⎠ ∂x where x, y, and z coordinates in the flow direction, across the film, and axial direction respectively (x, y, z in the equations refers to (xs,ys,zs) on fig 2.1a and b as the origins). Note all reference to), v1 and v2 are the runner and pad velocities in the transverse ydirection and F0 =  ∫  h  0  1  μ  dy ,  F1 =  ∫  h  0  y  μ  F2 = ∫  dy  h  0  F ⎞ y⎛ ⎜⎜ y − 1 ⎟⎟ dy μ⎝ F0 ⎠  (2.3)  Slider Bearings: Applying the following dimensionless group for slider bearings (the bar on top of each variable indicates dimensionless quantities as defined in the nomenclature)  u=  u ; U1  v=  y=  y ; h  z=  B v ; hi − ho U 1 z ; L  w=  h=  w ; U1  h ; hi − ho  x=  p=  15  x B  p (hi − ho ) μBU 1  2  (2.4)  and assuming U2 (stationary pad) and (v1-v2=0), the dimensionless form of the steadystate generalized Reynolds equation for slider bearings is [55] ∂ ⎧ 3 ∂p ⎫ 4 ∂ ⎧ 3 ∂p ⎫ ∂ ⎡ F1 ⎤ ⎨h F2 ⎬ + 2 ⎨h F2 ⎬ = ⎢h ⎥ ∂x ⎩ ∂x ⎭ λ s ∂z ⎩ ∂z ⎭ ∂x ⎣ F0 ⎦  (2.5)  where λ s = L B is the aspect ratio for sliders and F0 =  1  ∫μ 1  0  F1 =  dy ;  y  ∫μ 1  0  F ⎞ y⎛ ⎜⎜ y − 1 ⎟⎟dy 0 μ⎝ F0 ⎠  F2 = ∫  1  dy ;  yp  zp  Stationary pad  xp  U2=0  v(x,y)  v2  hi ys  (2.6)  u(x,y) h v1  zs  ho U1  xs  B  Runner  (a)  ϕ  Oil supply hole Oil groove  xb R1  xs  ys  W e  θ  R3  yb  R2 h  ω z x  fxs y  fys (b) Figure 2.1(a) Slider bearing (b) Journal bearing and notation  16  and pressure boundary condition p = p s at x = 0 ; p = 0 at z = 0 ; and ∂p ∂z = 0 z = ±1 / 2 at  (2.7)  Journal Bearings: Dimensionless parameters for journal bearings are x=  u =  x ; 2πR  u ; U  v=  y= R v ; CU  y ; h  z=  w=  w ; U  z L  h =  h ; C  P=  C 2P μi R 2 N s  (2.8)  These are substituted into the generalized Reynolds equation to give [55]  ∂ ⎧ 3 ∂p ⎫ 1 ∂ ⎧ 3 ∂p ⎫ ∂ ⎡ ⎛ F1 ⎞⎤ ⎨h F2 ⎬ + 2 ⎨h F2 ⎬ = ⎢ h ⎜1 − ⎟ ⎥ ∂x ⎩ ∂x ⎭ λ ∂z ⎩ ∂z ⎭ ∂x ⎣ ⎜⎝ F0 ⎟⎠⎦  (2.9)  where λ is the aspect ratio L/D and F0 , F1 and F2 are defined in equation 2.6, In addition to pressure boundary conditions equation 2.7, the Swifter-Stieber pressure boundary condition ∂p  ∂x  = 0 and p = 0 at x = xcav applies at the rupture zone.  Similarities in equations 2.5 and 2.8 for slider and journal bearings are clearly evident. Only minor adjustments are required to specify the length in the flow direction as x=2πR in journal bearings. The flow boundary condition on surface S1 (mid-plane) is zero which can be expressed mathematically as [56]  h 3 ∂p h − u ⋅n = 0 12 μ ∂x i 2  (2.10)  where xi takes the value of i=1,2 and corresponds to the coordinates, u is the velocity vector, and n is unit normal vector. It must be emphasized, however, that some of the gradient terms  ⁄  are not zero at the mid-plane as noted in section 4.3.2.  17  2.2.2 Energy and Viscosity Variation The three dimensional steady state energy equation is used to account for variations in film temperature.  ρc f u ⋅ ∇T − ∇ (k f ∇T ) = ρδTu ⋅ ∇p + f s  (2.11)  where the first term is convective, the second determines conduction, u is the velocity vector, kf the fluid thermal conductivity, cf the fluid heat capacity, ρ the fluid density, δ the compressibility, and fs represents source terms. Simplifying assumptions made are [2] 1. Conduction in the sliding (x) direction is negligible 2. Thermal conductivity and specific heat are constant, and 3. Velocity gradients in all but the transverse (y) direction are negligible. The form of the equation used in the finite element formulation is [2]  ⎡⎛ ∂u ⎞ 2 ⎛ ∂ 2T ⎞ ⎛ ∂T ∂p ⎞ ∂T ⎞ ∂T ⎛ ∂p ρc f ⎜⎜ u + w ⎟⎟ = k f ⎜⎜ 2 ⎟⎟ + ρδT ⎜ u + w ⎟ + μ ⎢⎜⎜ ⎟⎟ +v ∂z ⎠ ∂z ⎠ ∂y ⎝ ∂x ⎢⎣⎝ ∂y ⎠ ⎝ ∂x ⎝ ∂y ⎠  ⎛ ∂w ⎞ + ⎜⎜ ⎟⎟ ⎝ ∂y ⎠  2  ⎤ ⎥ ⎥⎦  (2.12)  in which the energy is expressed per unit volume. Dowson [54] and Ezzat and Rohde [57] show how to compute the velocity profiles u, v and w from the momentum and continuity equations (reproduced in Appendix A). (x, y, z in the equations refers to (xs,ys,zs) on fig 2.1a and b as the origins). Two methods of non-dimensionalising equation 2.12 are presented. The dimensionless parameter group the Pecklet (Pe), Prandtl (Pr), and Eckert (Ec) numbers present one method of making the energy equation non-dimensional  Ec =  U ; c f Ti  Pe =  ρc f Uho2 kf L  18  ;  Pr E c =  μ iU 2 k f Ti  (2.13)  where the subscript i indicates an input’s value at the inlet. These are substituted into equation 2.12, with the resulting equation ⎛ ∂T + ∂T + ∂T + +v +w ρ ⎜⎜ u ∂y ∂z ⎝ ∂x  ⎞ 1 ⎛ ∂ 2T + ⎞ Pr E c ⎛ ∂p ∂p ⎞ ⎟⎟ = ⎜⎜ ⎟+ ⎜⎜ u ⎟ +v 2 ⎟ ∂y ⎟⎠ Pe ⎝ ∂x ⎠ Pe ⎝ ∂y ⎠ 2 2 Pr E c ⎡⎛ ∂u ⎞ ⎛ ∂w ⎞ ⎤ + μ ⎢⎜ ⎟ + ⎜ ⎟ ⎥ Pe ⎣⎢⎝ ∂z ⎠ ⎝ ∂z ⎠ ⎦⎥  (2.14a)  where T+=T/Ti indicates the extent of temperature increase from inlet conditions. This formulation is applicable to both slider and journal bearing with little modification. However Jang et al. [55] derive a closed form of the energy equation  μ ⎛ ∂u ⎞ ∂T 1 ⎛ dh ⎞ ∂T κ1 1 ∂ 2T u + ⎜⎜ v − u y ⎟⎟ = 2 2 2 + κ1 2 ⎜⎜ ⎟⎟ ∂x h ⎝ dx ⎠ ∂y κ 2 h ∂y h ⎝ ∂y ⎠  2  (2.14b)  where two key parameters κ1 and κ2 referred to as temperature rise parameters appear and defined for slider and journal bearings as κ1 =  αμ i β U 2 ⎛ kfB  B ⎜ ⎜h −h o ⎝ i  ⎞ ⎟ ⎟ ⎠  2  , κ2 =U2  μi β  κ1 =  kf  αμ i β U 2 ⎛ R ⎞ kf R  ⎜ ⎟ ⎝C ⎠  2  ,κ2 =U2  μi β kf  (2.14c)  All terms can be determined from inlet conditions. The first temperature rise parameter κ1 is associated with viscous dissipation while κ2 relates oil properties and velocity. Slider Bearings: The heat flux at the oil film/pad interface is assumed equal kf  ∂T ∂y  = kp y b =1  ∂T p ∂y p  (2.15a) y p =0  while heat conduction in the pad is determined from ∂ 2T p ∂y p 2  +  ∂ 2T p ∂z p 2  =0  (2.15b)  19  where kf and kp are the fluid and pad thermal conductivities,  ,  and  the non-  dimensional coordinates of the pad as defined in the nomenclature, and T p = β (T p − Ti ) is the non-dimensional pad temperature, and β the temperature-viscosity coefficient. At the pad/ambient interface ∂T ∂y p  where  ,  y p =1  ⎛ = Bip ⎜ T p ⎝  y p =1  ⎞ − To ⎟ ⎠  (2.16a)  are the pad surface and ambient temperatures respectively, and Bip, the biot  number for the pad. At the film/runner interface heat flux is assumed equal kf  ∂T ∂y  = ks ys =0  ∂T s ∂y s  ys =0  (2.16b)  and the temperature determined from ∂ 2Ts ∂y s 2  +  ∂ 2Ts ∂z s 2  =0  (2.16c)  Journal Bearings: The energy equation as presented applies to journal bearings with modifications to account for heat transfer in the rupture zone. The rupture zone presents challenges as it comprises regions of alternating fluid and gas forming finger-like striations. Boncompain et al. [15] use two different energy equations in the rupture zone, one for the fluid and another for the gas sections. Critical to any modeling attempt is the identification of the rupture and reformation interfaces. It is expected that the shear rate decreases whereas the heat convection increases in the rupture zone while the lubricant is mixed with gas thus decreasing viscous dissipation significantly. Another approach to account for thermal characteristics of lubricant in this zone is to use the ‘effective length’ [15,58] which is defined as L’(x)=[hcav/h(x)]L and in dimensionless form as L 20  1 ⎤ ⎢h (x cav ) u (x , y , z )dy ⎥ dz 0 ⎣ 0 ⎦ L (x ) = 1/ 2 ⎡ 1 ⎤ ⎢h (x ) u (x , y , z )dy ⎥ dz 0 ⎣ 0 ⎦ 1/ 2 ⎡  ∫  ∫  ∫  ∫  for xcav < x < 1  (2.17)  where h(xcav) is the film thickness at the rupture interface and the other symbols have their usual meanings. Effective values of fluid conductivity k 'f , specific heat c 'f or ' viscosity μ are determined using a linear relationship of the form  Ψ’=ΨG+(ΨF-ΨG) L ( x )  (2.18)  where Ψ represents the property and the subscript denote gas (G) or fluid (F). Fluid film properties are assumed constant in the continuous sections in the transverse (y) direction at each position along the bearing length. The heat flux at the oil film/bush interface is kf  ∂T ∂y  = kb ys =1  ∂Tb ∂yb  (2.19a) yb =0  where kf and kb are the fluid and bush thermal conductivities and Tb = β (Tb − Ti ) the nondimensional bush temperature. Heat conduction in the bush is ∂ 2Tb ∂y b 2  +  ∂ 2Tb ∂z b 2  =0  (2.19b)  Convection applies at the outer surface of the bush and shaft. At the shaft end [15,75] ∂ Ts ∂z s  =− z s = ±1 / 2  Hs Bis ⎛⎜ Ts ⎝ ks  z s = ±1 / 2  − To ⎞⎟ ⎠  (2.20)  where Hs is the convective heat transfer coefficient at the shaft/ambient interface, Bis the Biot number for the shaft, and To the dimensionless ambient temperature. At the bush ends convection carries heat away from the shaft ends [15,75]  21  ∂Tb ∂z  =− z = ±1 / 2  (  Hb Bib Tb kb  z = ±1 / 2  − To  )  (2.21)  Inlet temperature conditions present some challenges as the fresh oil mixes with circulating oil in the bearing. The oil mixing model was verified inversely in that the oil temperature at the supply hole was assumed to be Ti and the model determines how the inlet conditions vary. A parabolic mixing model shown in Figure 2.2 for Dowson’s [59] data is determined, similar to the findings of other researchers [5,6]. Temperature variation across inlet film 41.00  Temp (C)  40.00 39.00 Temp  38.00 T = -5E+08z2 + 6358.5z + 312.4 R2 = 0.9983  37.00 36.00 0.0E+00  2.0E-05  4.0E-05 y-axis (m)  6.0E-05  8.0E-05  Figure 2.2 Inlet temperature variation for Dowson’s [59] bearing data. Computation of the inlet temperature was done following recommendations by Costa et al. [6] whereby the ratio of the oil groove width Lg to bearing width L determines the particular method used. If Lg/L>0.5 the energy conservation principle is applied at the oil groove whereas if Lg/L≤0.5 the energy conservation principle is applied to the section immediately downstream of the supply groove, over the full bearing width. Film pressure in the groove region for the case where oil groove is aligned with the load line as indicated in figure 2.1b was shown to be equal to the supply pressure. Simplifications made by Costa et al. [6] were used:  22  Tag=Ts+Cms(Tin-Ts)  for Lg/L>0.5  Tag=Tin  for Lg/L≤0.5  Tig=Tag  for all Lg/L  (2.22)  where Tag is the temperature of the oil coming out of the groove in the axial direction, Tig is the temperature of the oil in the bearing clearance on the groove region, and Cms is an empirical mixing coefficient between zero and one (0.5 was used as recommended by Costa et al.). If film reformation was detected upstream of the groove, Ti =  Ti =  Tup Qup + Trv Qrv + [Qs − (1 − C ms )Qag ]Ts Qin + C ms Qag  Tup Qub + Ts Qs s  for a/b>0.5  for a/b≤0.5  Qub + Qs  (2.23)  If the film reformation was over the groove region or downstream, the oil inlet temperature was determined by Ti =  Tup Qup + Trv Qrv + [Qs − (1 − C ms )Qag ]Ts Qin + C ms Qag  Ti =  Tup Qre + Ts Qs  for a/b>0.5  for a/b≤0.5  Qre + Qs  (2.24)  In equations 2.23 and 2.24 Tup is the average lubricant temperature of the re-circulating oil flow, Trv is the average lubricant temperature of the reverse flow, Qs is the oil supplied to the bearing, Qre is the re-circulating flow coming from the rupture region, Qup is the oil flow calculated upstream of the groove for the groove axial extent, Qub is the oil flow upstream of the groove calculated for the full bearing width, Qin the flow rate at the inlet section, and Qrv is the reverse flow.  23  Viscosity and temperature equations form the final couple of equations in the set. They are related to temperature variations by  μ = μi exp[− β (T − Ti )]  (2.25a)  ρ = ρ o [1 − δ (T − Ti )]  (2.25b)  where β is the temperature-viscosity coefficient and μi the inlet viscosity and δ the temperature-density coefficient. The viscosity equation couples the Reynolds and energy equations and can be expressed in dimensionless form by introducing the dimensionless temperature T = β (T − Ti ) which, when substituted into 2.25 gives  μ = exp[− T ]  (2.26a)  ρ = 1 − δ (T − 1)  (2.26b)  Equation 2.26a is based on the effective viscosity concept [55] (as defined in the List of Symbols) while 3.26b eminates from equation 2.14a δ =  Pr Ec δTi Pe  (very small). The  viscosity-temperature relationship has a limited range of applicability - temperature variations on the order of a few tens of degrees – thus will not be applied to peak temperature predictions well in excess of 100°C.  The last two sections present literature review on methods of solving the thermohydrodynamic problem in journal bearings. Critical in this endeavor is proper specification of boundary conditions. Although not always easy to do, specification of boundary conditions at the groove/oil inlet is critical.  2.3  Finite Element Formulation with Streamline Upwind Petrov-Galerkin  The slider in figure 2.3 is shown for reference (only half the bearing is used due to  24  symmetry). Mid-plane surface, S1 is shown hatched. Numerical solution of the pressure problem is based on the variational/weak formulation. Assume the domain of interest (shown in figure 2.4) is bounded in R3 with a piecewise boundary Γ and let n( x , y, z ) denote the spatial coordinates of the outward unit vector of a generic point on Ω . Let a portion of the boundary be Γ1 on which pressure computations are performed. The pressure p is the solution of the minimization problem [60] 2 2 h 3 1 ⎡⎛ ∂p ⎞ ⎛ ∂p ⎞ ⎤ ⎛ ∂p ∂p ⎞ min Π (p ), with Π ( p ) = ∫ ⎢⎜ ⎟ + ⎜ ⎟ ⎥ − h ⎜ + ⎟dΓ p ⎝ ∂x ∂z ⎠ ⎢⎝ ∂x ⎠ ⎝ ∂z ⎠ ⎦⎥ Γ1 6 μ 2 ⎣  (2.27)  As mentioned earlier the Galerkin weighted residual formulation is not reliable in convective dominated flows as it leads to spurious node-to-node oscillations in the solution. Some form of upwinding is needed, and the streamline upwind Petrov-Galerkin (SUPG) shows good stability in initial and boundary value problems and can handle  ~ and backflow at the inlet [61]. The energy equation is multiplied by a test function w integrated over the domain Ω to get the weak formulation. ⎡ ⎛ ∂T ∂T ∂T ⎞ 1 ⎟⎟ − +v +w ⎢ ρ ⎜⎜ u ∂ ∂ ∂ x y z ⎝ ⎠ Pe ~⎢ w 2 2 ⎢ ∫Ω ⎡ ⎤ ⎢ − Pr E c μ ⎢⎛⎜ ∂u ⎞⎟ + ⎛⎜ ∂w ⎞⎟ ⎥ ⎜ ⎟ ⎜ ⎟ ⎢ Pe ⎢⎣⎝ ∂y ⎠ ⎝ ∂y ⎠ ⎥⎦ ⎣  The second term in equation 2.29,  ⎛ ∂ 2T ⎜⎜ 2 ⎝ ∂y  ∫ w(1 P )(∂T ~  e  ⎞ Pr E c ⎛ ∂p ∂p ⎞ ⎤ ⎟⎟ − +w ⎜u ⎟⎥ ∂z ⎠ ⎥ Pe ⎝ ∂x ⎠ (2.29) ⎥dΩ = 0 ⎥ ⎥ ⎦ 2  ∂y 2 )dΩ is integrated by parts to reduce  Ω  ~ consists of two parts: the classical the derivatives by one order. The weight function w Galerkin weighting function w and a discontinuous part ζ  (the streamline upwind  ~ becomes w ~ = w + ζ . The parameter) which is added so that the weighting function w  25  finite element weak formulation is solved by approximating the solution by a finite number of basis functions and approximating the weighting function w by the same set. L/2 Inlet  B  yp zp xp  S3 A  hi  hp  ys zs  S1  xs  B  U  ho  Outlet S5  Figure 2.3 Three dimensional slider bearing geometry. Surface S1 is the mid-plane, S2 and S5 are the inlet and outlet surfaces. n( x , y, z )  Γ R3  Ω R2 R1  Γ1  (a)  Projection of nodal pressure values on Γ1 into volume domain Ω  Pressure calculated on bottom surface  (b) Figure 2.4(a) Solution domain Ω showing boundary Γ1 and unit normal vector n. (b) Domain Ω transformation (only a small part of the transformation is shown).  26  This results in the standard finite element interpolation or shape functions Ni which doesn’t yet account for the extra part in the weighting function. The ζ-part of the basis function is defined at the element level as ζe [62]  ζe =  Δxξ u ⋅ ∇N i 2 u  (2.30)  where Δx is the size of the element in the direction of flow, ξ is a parameter defining the type of upwinding (in the classical upwind scheme ξ=1), ∇Ni is the gradient of the interpolation (shape) functions, and ||u|| is the square root of the inner product of u, (||u||=(uiui)1/2). In Appendix E a comparison of thermal prediction results from Galerkin and SUPG is explored. In the finite element method approach the lubrication problem is discretized into a finite number of unknowns. This is achieved by dividing the solution domain Ω into hexahedral (brick) elements and expressing the field variables (pressure, velocity and temperature) in terms of approximating trial functions (Ni) within each element. The volumetric domain  Ω is divided into Nve linear hexahedral elements. The field variables (velocity and temperature) are expressed in terms of approximating trial or shape functions (Ni) within each element whereas the boundary surface Γ1 is discretized into Ns bilinear rectangular elements such that [56] Ω  Ω,  Ω  (2.31a)  Γ  Γ,  Γ  (2.31b)  where φ1 represents the volumetric solution field and φ2 the surface solution vector, and and  represent the union and intersect of the elements. The discretized field  variables are given by [52,53] 27  u ≈ e  N ne  N ne  N ne  N se  k =1  k =1  k =1  k =1  ∑ N kv u k , v e ≈ ∑ N kv vk , Te ≈ ∑ N Tk Tk , and p e ≈ ∑ N kp pk  (2.32)  where Nne=8 for the brick elements and Nse=4 for the bilinear rectangular elements and ,  and  are the shape functions for velocity, temperature and pressure  respectively. These elemental discretizations are introduced into equation 2.27 for the pressure [52,53] ⎧ 3 ⎪⎪ h ⎨ Γ1e ⎪12μ ⎪⎩  ⎡ ⎛ ⎢∂ ⎜ ⎢ ⎜ ⎢ ∂x ⎜⎝ ⎣  ∫  N se  ∑ j =1  2  ⎞ ⎛ ∂ ⎜ ⎟ p ej N jp ⎟ + ⎜ ⎟ ∂y ⎜ ⎠ ⎝  N se  ∑ j =1  ⎞ ⎟ p ej N jp ⎟ ⎟ ⎠  2⎤  ⎛ ⎥ ∂ ⎜ ⎥−h ⎜ ∂x ⎜ ⎥ ⎝ ⎦  N se  ∑ j =1  ⎞⎫⎪ ⎟⎪ p ej N jp ⎟⎬dΓ1 = 0 ⎟⎪ ⎠⎪⎭  (2.33)  Equation 2.33 is minimized with respect to the elemental field variable p ej and on further simplification becomes [52,53] N se  ⎡h 3 ⎢ Γ1e ⎢ 6μ ⎣  ∑∫ j =1  ⎛ ∂N jp ∂N p ∂N jp ∂N p ⎞⎤ i ⎟⎥ e i ⎜ + pj − ⎜ ∂x ∂x ∂z ∂z ⎟⎥ ⎝ ⎠⎦  ∫  Γ1e  h  ∂N ip dx = 0 ∂x  (2.34)  The energy equation can be discretized in a similar manner to yield [52,53,63] ⎧ ⎛ ⎞ v ⎛ N ne ⎛ N ne ⎞ v ⎛ N ne ⎞ ⎞ ⎡⎛ N ne ⎞ v ⎤⎫ ⎪~ ⎜ e v ⎟ ∂N j ⎥ ⎪ e v ⎟ ∂N j v e e v ⎟ ∂N j ⎜ ⎜ ⎟ ⎟ ⎢⎜ ⎜ + − − 1 1 β + w N T u N v N w N ⎨ i⎜ ⎬ k⎟ k⎟ k k k⎟ k ⎜ ⎟ ⎟⎟ ⎢⎜ ∂x ⎜ k =1 k ∂y ⎜ k =1 k ∂z ⎥ ⎪ Ωe ⎪ ⎜ ⎠ ⎝ ⎝ k =1 ⎠ ⎠ ⎠ ⎣⎝ k =1 ⎝ ⎠ ⎦⎭ ⎩ ⎝  N ve  ∑∫ j =1  ∑  ∑  ∑  ∑  ⎛ ⎛ ⎛ N ne ⎞⎞⎞ v e ⎟⎟⎟ ~ Pr Ec ⎜ exp⎜ − ⎜ w N T i k k ⎜⎜ ⎜ ⎟ ⎟⎟ ⎟⎟ Pe ⎜⎜ Ωe ⎠⎠⎠ ⎝ ⎝ k =1 ⎝ 2 2 ⎡ ⎛ ⎛ N ne ⎞ ⎛ N ne ⎞ ⎞⎤ ⎢ ∂ ⎜⎜ e v⎟ e v ⎟ ⎟⎥ ⎜ + u N w N k⎟ k ⎟ ⎟⎥ ⎢ ∂y ⎜ ⎜ k k ⎜ ⎢ ⎜⎝ ⎝ k =1 ⎠ ⎝ k =1 ⎠ ⎟⎠⎥⎦ ⎣ ⎡⎛ N ne ⎞⎛ ⎛ N ne ⎞ ⎞ ⎛ N ne ⎞⎛ ⎛ N ne ⎞ ⎞⎤ ⎟ e v ⎟⎜ d ⎜ v e ⎟⎟ ⎜ e v ⎟⎜ d ⎜ ⎜ ⎢ + u k Nk ⎜ N k pk ⎟ + wk Nk ⎜ N kv pke ⎟ ⎟⎥ ⎟⎜ dx ⎜ ⎟⎟ ⎜ ⎟⎜ dz ⎜ ⎟ ⎟⎥ ⎢⎜ ⎠⎝ ⎝ k =1 ⎠ ⎠ ⎝ k =1 ⎠⎝ ⎝ k =1 ⎠ ⎠⎦⎥ ⎣⎢⎝ k =1  −  1 Pe  ⎡⎛ ∂W ∂N v ⎞ ⎤ j ⎟ ⎢⎜ i ⎥T je = ⎢⎜⎝ ∂y ∂y ⎟⎠ ⎥ ⎣ ⎦  ∑  ∑  ∑  ∫  ∑  ∑  ∑  ∑  (2.35) where Nne is the number of nodes for each elemental, N is the shape function of the kth v k  ~ are the Galerkin and SUPG weighting functions respectively volume element, wi and w i and equation 2.26 has been substituted for ρ and μ while the subscript e represents a variable at the elemental level. One important observation about the discretization in  28  equation 2.35 is that the gradient pressure vectors ∂p / ∂x and ∂p / ∂z on boundary Γ1 have to be extended into the volume domain Ω in order to compute the velocity components u and v (figure 2.4b). The governing equations are used to develop the model using the finite element program Sepran [62], an FE program developed at Delft University. The program is structured into ‘standard problems’ that mainly encompass solution to secondorder elliptic and parabolic differential equations. One advantage the program has is the ease with which user subroutines are implemented in the Fortran programming language. An example is the computation of compressibility terms in equation 2.35, velocity components u, v, and w, effective viscosity μeff and density ρ. Performance parameters are the load capacity, frictional drag force, and side flow. The non-dimensional load carrying capacity and frictional force for a slider are [9] 1 1⎡  ∂ ⎧ 3 ∂p ⎫ 4 ∂ ⎧ 3 ∂p ⎫ ∂ ⎡ F1 ⎤ ⎤ ⎢ ⎨h F2 ⎬ + 2 ⎨h F2 ⎬ − ⎢h ⎥ ⎥ dx dz ∂x ⎭ λ s ∂z ⎩ ∂z ⎭ ∂x ⎣ F0 ⎦ ⎥⎦ 0⎣ ⎢ ∂x ⎩  ∫∫ 0  (2.36a)  (2.36b) is determined in Appendix A and  where velocity component  is the  differential of equation A_2 with respect to y. Shaft frictional torque Tfr in journal bearings is determined by integrating shear stresses (  /  at the shaft and bush interfaces [54] (2.37)  where velocity components  and  are presented in Appendix A. The first term is due  to film pressure and second due to velocity-induced shearing. Load capacity for journal  29  bearings is presented in section 2.4. Equations 2.36 and 2.37 were solved using user subroutine which performs the integrations using current values of viscosity μ . The flow rate has a pressure and velocity-induced component [64] (2.38) where  is determined in appendix A. Side flow in slider bearings can also be determined  by assuming conservation of mass flow and finding the difference between the inflow and outflow, QL = Qi − Qo .  2.4  Modeling Procedure  The equations presented in sections 2.2 and 2.3 form the core of the bearing template developed. A flow chart of the implementation is shown in figure 2.5. Current values of viscosity μ are used to compute the pressure field from Reynolds equations (2.5 and 2.8) with the boundary conditions equations 2.7 and 2.9. A conjugate gradient iteration method is used. The pressure is used in the next step where the Reynolds equation is solved with the Swift-Stieber boundary condition using an over-relaxation technique with constraint p≥0. Two important forces on the bearing are the load capacity as determined from the Reynolds equation Wd and the actual load W. A force balance must exist between the applied load and load capacity. The procedure used to achieve this is a. Assume an initial value for eccentricity e, compute eccentricity ratio (ε=e/c) (normal range 0.5≤ε≤0.9) as well as attitude angle ϕ and compute the film thickness h(θ) from (see figure 2.1b) h(θ)=c(1+εcos(θ-ϕ))  30  (2.39a)  where c=(R1-R2) is the radial clearance, θ and ϕ the angle from the load line and attitude angle respectively. The minimum film thickness is located at θ=π+ϕ hmin=c(1-ε)  (2.39b)  b. Solve the Reynolds equation for pressure and load-carrying capacity Wd. c.  For this capacity compute the reaction load fR=(fx,fz)T and attitude angle from fx =  ∫  Ω  − p cos(ϕ )dΩ  (2.40a)  f z = ∫Ω − p sin (ϕ )dΩ  and  (2.40b)  fR  (2.40c)  and attitude angle ϕ which is ⎛ fy ⎞ ⎟ ⎝ fz ⎠  ϕ = tan −1 ⎜  (2.40d)  Sommerfeld number can be computed from  S=  μeff N s DL ⎛ R ⎞ W  ⎜ ⎟ ⎝C ⎠  2  (2.40e)  d. If fR<W increase the value of eccentricity ratio ε, if fR>W decrease its value. e. Compute the attitude angle from equation 2.40d. If ϕ<ϕa (where ϕa is value from step a) decrease value of ϕ (in step a), or if ϕ>ϕa increase value of ϕ (in step a). f. Go back to step a and repeat procedure. g. When convergence is obtained on both eccentricity ratio and attitude angle, the flowchart in figure 2.5 continues at stage 4 by computing velocity components u, v, w.  31  h. modify the attitude angle until the direction of the calculated load from 5c is equal to the imposed load direction. i. Once equilibrium is achieved the flowchart in figure 2.5 continues at stage 3 by computing pressure gradients dp/dx and dp/dz followed by velocity components u, v, w. The equations for the velocity fields are numerically integrated in a user subroutine using the trapezoidal rule. Velocity component w is determined from the continuity equation A-7 (Appendix A). Once the velocity field is determined, the energy equation is solved for temperature, which is used to compute new density and viscosity fields (note that density variation is so small it may be ignored). Computation of the thermal field in the rupture zone is achieved by considering effective length equation 2.17 and the linear relationship in 2.18 for other physical properties (thermal conductivity and heat capacity). New values of effective viscosity are determined with the thermal field and compared to previous values. The effective viscosity is the average viscosity value over the film nodes determined using equation 2.26a and dividing by the number of nodes. If this value is more than a predetermined value εr (set to μ≤0.01 or if the number of iterations iter exceeds a specified limit Lm (set to 100) in case of slow/no convergence) the iteration repeats (outer loop in figure 2.5) otherwise the program proceeds to computing bearing operating parameters  ,  ,  . ε r ≤ μ i − μ i −1  (2.42)  Jang and Khonsari [55] recommend an empirical equation for determination of ε that may be used to verify equation 2.40 (note that this is used for verification purposes only and used at the end of the process).  32  ε = −0.36667 +  0.49737 L/D  − 0.20986 ln (S )  for 0.5 ≤ L/D ≤ 1.5  (2.41)  Periodical boundary conditions are applied to surfaces S3 and S5 in figure 2.3. Similar conditions are set for the bush and shaft ends. A typical run using data from Dowson et al. [59] is given below a. Initial guesses: ε=0.5, ϕ=0.9 radians. b. Initial results from solving Reynolds: horizontal component of load fx: 2.81838E+03 vertical component of load fy: 4.94780E+03 modulus of load fR: 5.69421E+03 attitude angle: 1.05301E+00 c. W=11 000N, therefore W>fR increase ε to 0.64. Second iteration gave the results horizontal component of load 7.05888E+03 vertical component of load 8.51415E+03 modulus of load 1.10598E+04 attitude angle 8.78576E-01 Third iteration with ε=0.63: horizontal component of load 7.24169E+03 vertical component of load 8.28417E+03 modulus of load 1.10032E+04 attitude angle 8.52443E-01 The modulus of load or fR=11 003 ≈ W; the iterations can be stopped. The flowchart continues at stage 4 in figure 2.5.  2.5  Chapter summary  The streamline upwind Petrov-Galerkin finite element method is formulated for threedimensional thermo-hydrodynamic lubrication problems. Governing equations and boundary conditions in the critical areas of the groove and rupture zone are presented. The model is used to simulate slider bearings thermal behavior in chapter 3 and journal bearings in chapter 4.  33  Start. Input parameters Ti, μa,, n, k,ρ, D, C  Use current value of μ to compute p using constrained minimization and over-relaxation. Determine ε, h, ϕ (eqns. 2.39-2.40) at force balance equilibrium  Find gradients ∂p/∂x, ∂p/∂y Compute the velocity components u, v, w using Navier-Stokes and continuity equations Solve energy equation iteratively for temperature field T in film. Use conduction eqn to determine temperature in pad/bush  Compute new viscosity and density fields  Compute new pressure field pnew with updated viscosity , frictional torque Tr, power loss PL,, QL for current values of viscosity  No  Is |µnew-µold|≤εr Or is Iter>Lm  Yes Compute output parameters Tr , QL and solve for conduction in pad/bush  Yes  Change value of θ?  No End  Figure 2.5 Flowchart for thermo-hydrodynamic lubrication model.  34  Chapter 3 Application of SUPG FE to Slider Bearings  3.1 Introduction The streamline upwind Petrov-Galerkin (SUPG) finite element method formulated in chapter 2 is chosen for its robust nature and rapid convergence in convection dominated flows. The method is applied to slider bearings analysis and design in this chapter. Slider bearings are widely used in applications such as mechanical seals, collar thrust bearings, machine tool guides, and piston rings. They have good load-carrying capacity, excellent stability and durability. A tilted pad slider bearing consists of two surfaces as shown in figure 2.1a separated by a lubricant film. The direction of the runner motion and inclination of the planes are such that a converging oil film is formed between the two surfaces and the physical wedge pressure-generating mechanism is developed in the lubricant film. This enables the bearing to support a load. The flow is viscous, leading to heat generation, dissipation and increase in film temperature.  3.2  Lebeck Slider Bearing Data  The test case presented by Lebeck [65] and used by Kumar [9] in a two-dimensional model with input values in Table 3.1 is used. Mesh density sensitivity analysis was performed with the thermo-hydrodynamic SUPG model; results appear in Table 3.2. The maximum number of pressure and temperature iterations for convergence is also indicated. An analytical load capacity is included and used as a benchmark. A 30x30x10 (x-y-z) mesh is finally chosen. The length (x) and height (y) values are similar to those used by Kumar [9] while an axial value of ten elements is a good compromise in  35  computation time and accuracy. The mesh size was chosen after trials with other densities – a compromise had to be made on mesh density versus computing time.  Table 3.1 Input values for model [65]. Input Parameter  Value  Sliding Speed, U Leading edge height, hi Bearing Length L Pad thickness hp Lubricant heat capacity cf Lubricant conductivity kf Pad conductivity kb Runner conductivity ks Inlet temperature Ti Inlet lubricant density ρi Inlet viscosity μi Viscosity-temp coeff. β Density-temperature coeff. δ Ratio k (ho/hi) Inlet pressure pi (gauge)  20 m/s 5x10-5 m 0.1 m 0.1m 1926 J/kg K 0.132 W/m K 60 W/m K 50 W/mK 310 K 897.1 kg/m3 0.0174 Pa.s 0.035 /K 0.0012 /K 0.4 0.0  Table 3.2 Mesh sensitivity analysis Mesh size Load (x,y,z) Capacity 5x5x5 0.095 10x10x10 0.106 20x20x10 0.108 0.4 30x30x10 0.109 40x40x10 0.109 50x50x10 0.109 Analytical: 0.109 Max. no. of pressure iterations Temperature iterations k  Diff. (from analytical model) % 12.74 2.94 0.66 0.22 0.07 0 6 10  3.2.1 Model Convergence Performance The iterative procedure in figure 2.5 showed convergence for load capacity computations of less than ten iterations for L/B=10 and less than five for L/B=1 as shown in figure 3.1.  36  Table 3.3 has the limiting values for the outer loop in figure 2.5 and an indication of the convergence values. The power and advantage of SUPG is evident in the finite bearing case where, at L/B=1 convergence takes place after only three iterations. A similar pattern is evident in journal bearing simulations.  Load Capacity as Function of Number of Iterations 1.40E+06  1.20E+06 L/B=10 L/B=5 L/B=3 L/B=2 L/B=1  Load Capacity, N/m  2  1.00E+06  8.00E+05  6.00E+05  4.00E+05  2.00E+05  0.00E+00 1  2  3  4  5  6  7  8  Number of iterations  Figure 3.1 Convergence of load capacity solution.  Table 3.3 Convergence criteria Maximum number of iterations – outer loop Percentage variation – 1st to 2nd iteration Percentage variation – 5th to 6th iteration Pressure iteration: max no. Temperature iterations:  37  15 57.9 4.9 6 10  9  10  3.2.2  Predicted Pressure Values  The effect of the width to length ratio L/B on predicted pressures generally indicates whether a bearing is infinite or finite. Initially adiabatic conditions on upper and lower surfaces and an inlet temperature Ti=To are assumed. The 3-D SUPG model formulated in section 2.5 is used to predict peak pressure and load capacity. Results appear in table 3.4 and figure 3.2. Values of peak pressure for L/B ratios of ∞ (2-D model), 20, 10, and 5 are coincident as expected [66]. Percent difference in table 3.4 indicates model predictions for the infinitely wide (2-D) and 3-D sliders as a function of L/B ratio for both peak pressure and load capacity. Both fall dramatically at L/B values less than 2. What is significant here is the sudden drop in load carrying capacity from L/B of 2 to 1 (square slider). A 3-D model should be used for L/B less than 3. Table 3.4 Isothermal peak pressure and load capacity, k=0.4 L/B ratio 2D Model 20 10 5 2 1 0.5  Peak Pressure (p) 0.2564 0.2564 0.2564 0.2562 0.2378 0.1669 0.0806  Difference (%) 0.00 0.00 0.09 7.3 34.9 68.6  38  Load capacity  W 1.572 1.572 1.572 1.571 1.453 0.9996 0.4437  Difference to 2D model (%) 0 0 0.06 7.5 36.4 71.8  SUPG Predictions as function of L/B ratio 0.30  p_2D  0.25  p_20:1 p_10:1 p_5:1  0.20  p_2:1 p_1:1 p_0.5:1  0.15  0.10  0.05  0.00 0.00  0.10  0.20  0.30  0.40  0.50  0.60  0.70  0.80  0.90  1.00  x/L  Figure 3.2 Non-dimensional isothermal pressure predictions as a function of L/B. The first four lines (L/B=∞, 20, 10, 5) coincide. Arnell [66] suggests sliders are considered finite when L/B<2 which corresponds to a peak pressure difference of 7%. Of more significance is the thermo-hydrodynamic effect on peak pressure predictions as shown in Table 3.5 which indicates the difference between isothermal and THD values for several L/B ratios. Consideration of the thermohydrodynamic effect must be made at the bearing speed tested (20 m/s).  39  Table 3.5 THD peak pressure and load capacity comparisons L/B ratio 2D Model 20 10 5 2 1 0.5  Peak Pressure Isothermal (p) 0.2564 0.2564 0.2564 0.2562 0.2378 0.1669 0.0806  Peak Pressure THD (p) 0.144 0.144 0.144 0.144 0.135 0.096 0.045  Difference 0.112 0.112 0.112 0.112 0.103 0.071 0.036  Load Capacity, Isothermal  Load Capacity, THD  W  W  1.572 1.572 1.572 1.571 1.453 0.9996 0.4437  0.950 0.895 0.439 0.208 0.068 0.022 0.005  Difference 0.622 0.677 1.13 1.36 1.39 0.978 0.439  3.2.3 Thermal Profile Thermal fields in the film from the two-dimensional [9] and SUPG models for L/B=10 are compared in figure 3.3. The isotherms follow a similar pattern in both for this aspect ratio. The scale of the y axis representing the film thickness is magnified 1000 times since it is much smaller than the other two dimensions (length L and width B). Adiabatic conditions and an inlet temperature Ti=310°K and L/B=10 are used similar to Kumar’s model for comparison. Peak temperature prediction for the two dimensional model is 75.4°C (normalized temperature T+=1.124) at the upper corner of the trailing edge while that for the three dimensional model is 74.7°C (T+=1.121) at the same location. These results assume adiabatic conditions. Figure 3.4 is a 3D model of the film showing the predicted thermal field. Maximum temperatures occur at the trailing edge on the mid-plane surface (S1 in figure 2.3,  X = 1, Y = 1, Z = 0 ) at the film/pad interface. Peak temperature reduces to 66.3°C for L/B=0.5 indicating the increased effect of heat carried out through leakage flow (figure 3.5). Peak temperature is at the mid-plane surface. Table 3.6 shows that the fraction of leakage flow increases from 7 to 50% as L/B reduces from 10 to 0.5. 40  More industrially applicable conditions, such as those used by Kumar [9]: p i = 1.0 , Ti + = 1.2(372 K ) , T p+ = 1.4(434K ) and Ts+ = 1.0(310K ) are used to run the model. These  normalized temperatures are relative to the inlet temperature in °K, all other variables are in Table 3.1). Inflow boundary conditions were checked and updated as well as using a finer grid mesh N=30x30x30 with improved results as shown in figure 3.6a. A clearer picture emerges when predictions are compared to Kumar’s results by digitizing and superimposing these onto figure 3.6a. Figure 3.6b shows a couple of isotherms. What is evident from this is that although relatively close, it is difficult to reproduce Kumar’s results exactly. Isotherms generally follow similar patterns in both cases although model predictions appear to be shifted slightly upwards especially at the lower regions.  41  (a)  (b) Figure 3.3(a) Normalized temperatures for 2-D and (b) 3-D SUPG models showing general patterns of the isotherms.  42  1 0.8 0.6 0.4 0  0.2  Y+  Figure 3.4 Contour map of the thermal field, L/B=10 (hi=50μm). The arrow shows flow direction.  43  1 0.8 0.6  44  Y+  0.4 0.2 0  Figure 3.5 Thermal field with L/B=0.5 (hi=50μm).  (a)  (b) Figure 3.6(a) Model isotherms for Lebeck’s data [65] (b) Comparison of temperature profiles from SUPG model (solid lines) and Kumar [9] (dashed lines).  45  3.2.4  Leakage Flow Pattern – Flow Gradient Method  The SUPG model is used to predict flow patterns as a function of L/B in more detail. Table 3.6 shows the average side flow as a function of L/B, also shown in figure 3.7 for 0.1<k<0.9 and different L/B ratios. Flow is determined at k intervals of 0.1 and averaged out. Boxed numbers represent the ratio of side leakage to inflow as a percentage. It increases significantly for L/B ratios less than two. The values presented in table 3.6 and figure 3.7 are for average dimensionless flow rates for the range 0.1<k<0.9 and represents the general behavior of slider bearings as a function of aspect ratio L/B. Figure 3.8 follows from the general observations evident from figure 3.7 and is applicable to slider bearings in general.  The gradient in figure 3.8 is determined manually at each point on figure 3.7. A differencing technique is used to compute the change in tilt ratio d(L/B) and flow  .  The gradient gives an indication of how the flow is changing as a function of length to width ratio L/B. The minimum point IP represents, for a given tilt angle k, rate of change of leakage flow as a function of L/B ratio is zero. This point is at L/B=2 in figure 3.8. Lower values of the L/B give more leakage flow but far less load-carrying capacity whereas higher L/B values give more load capacity but less leakage flow and consequently higher operating temperatures. Thus Ip gives a good compromise between load capacity and leakage flow for optimum operating temperature. As pointed out earlier higher leakage flows are beneficial in lowering peak operating temperatures. The gradient of side flow ∂QL/∂(L/B) is high from zero to L/B=1, decreases rapidly between L/B=1 and L/B=2 to almost zero, gradually increases from two to ten, then becomes nearly constant above ten. It decreases rapidly from zero to one. It can be argued that the influence of side flow is most significant between zero and two. Even though load capacities are relatively low in this region, the beneficial effects of lower 46  operating temperatures are realized. Another interesting observation is that even at L/B=60 the percentage of side flow to inlet flow is 4.9 although the slider is in the infinite range. Figure 3.8 illustrates the Flow Gradient Method of determining optimum leakage flow by computing the gradient dQ / d (L / B) for 0.1<k<0.9. The graph bottoms out at a ratio of two, increases up to a ratio of ten and remains nearly constant thereafter. To the best of our knowledge this work represents the first definition of the Flow Gradient Method.  Table 3.6 Dimensionless flow rate as a function of L/B1  L/B 60 40 20 10 5 2 1 0.5 0.25  ,  , and  Qv 46.4 30.9 15.5 7.86 4.05 1.79 1.02 0.57 0.30  Qo  QL  44.1 29.4 14.7 7.29 3.60 1.37 0.63 0.29 0.13  2.3 1.6 0.87 0.57 0.46 0.42 0.39 0.28 0.17  Ratio of side leakage to inflow, % 4.9 5.1 5.6 7.2 11.2 23.3 38.0 49.7 55.4  are the inflow, outflow and side leakage flow respectively  47  Dimens sionless leakkage flow as function of L/B, L 01<k<0..9 2..5 4.9 2..0  1..5  5.1 -  1..0  0..5 3 38.0  5.6  23.3 11.2 7.2 49.7 55.4  0..0 0  10  20  30 40 gth to width ratio L/B Leng  50  60  Figure 3.7 Diimensionlesss side leakagge Q L as a fuunction of L//B, 0.1<k<0.9. Boxed numbeers represent the ratio of side leakagee to inflow as a a percentaage.  Flow Gradient G a Functio as on of L/B,, 0.1<k<0.9 0.5 0.45 0.4 0.35 0.3 0.25 dQ/d(L/B)  0.2 0.15 0.1 0.05 0  0 2 4 6 8 10 12 14 1 16 18 20 0 22 24 26 28 2 30 32 34 4 36 38 40 L/B  Figgure 3.8 Flow w gradient, 0.1<k<0.9.  48  In Figure 3.9 the non-dimensional side leakage pattern from point A to B in figure 2.3a is shown as a function of L/B for k=0.4. Numbers in boxes indicate the width to length ratio L/B. Side leakage increases from L/B=20 to a peak at L/B=1 and L/B=2 then declines as L/B decreases. This is an important observation, as it implies the bearing is optimized in  terms of side leakage in the range L/B=1-2.  Side leakage profiles 0.25  2 1 3  0.20 Non-dimensional leakage flow  L/B=20 L/B=10 L/B=1 0.5  L/B=0.5  0.15  L/B=2 L/B=3 10  0.10  L/B=20  0.05  0.00 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  Figure 3.9 Side leakage as a function of L/B, k=0.4.  Although the fact that finite sliders run at significantly lower temperatures than higher width ones has been long established [67] the effect on the point of the most intense viscous dissipation has not been widely reported. Peak temperature increase lowers lubricant viscosity which leads to lower load carrying capacity. Slider geometry can be  49  adjusted or new geometries can be tried in order to reduce the effect of higher temperatures on load carrying capacity. The effect of tilt ratio k on peak temperatures and load carrying capacity are shown in figure 3.10 for L/B of 20, 10, 5, and 1. All variables not specified are in Table 3.1. Boxed numbers show peak temperature predictions at each k value. A high tilt angle (low k value) results in high peak temperatures. The highest load capacity occurs at k=0.5 in both cases. Even though high values of k have a lower peak temperature, load capacity is less. A tilt ratio of 0.5 seems ideal for slider design, although slightly higher values (k=0.6-0.7) are preferred for their lower peak temperatures. The high peak temperatures at L/B=10 and 20 and tilt ratio k=0.1 suggest high risk of metal to metal contact for these values.  Load Capacity as a Function of k 1.00 69 C 64 C  0.90  L/B=10 L/B=1  75 C  L/B=5  61 C  0.80  L/B=20  0.70 82 C  0.60  58 C  NonDim Load Capacity  0.50 69 C  64 C  75 C  0.40  61 C 56 C  82 C  0.30  58 C  0.20  69 C  75 C  55 C  64 C  83 C  0.10  56 C  61 C 58 C  56 C 87 C  77 C  71 C  66 C  0.2  0.3  0.4 0.5 Parameter k  62 C  59 C  57C  56 C  49 C  0.00 0  0.1  0.6  0.7  0.8  0.9  1  Figure 3.10 Load capacity as a function of k. Non-dimensional load capacity is as defined in nomenclature.  50  3.3  Chapter Summary  The three dimensional SUPG thermo-hydrodynamic model developed is used to analyze slider bearing characteristics at different width to length and tilt ratios. Pressure predictions for the non-isothermal case are significantly less than the isothermal case as expected. Model predictions show square sliders running at lower temperatures than wider ones. This is due to the higher side leakage; however their load carrying capacity is significantly less. The significance of side leakage on peak temperature is evident. The model also shows that side leakage peaks in the region L/B=1-2. The Flow Gradient Method to predict optimum flow conditions is proposed. Even though high values of the tilt ratio have a lower peak temperature, the load capacity is less and an tilt ratio of 0.50.7 seems ideal for slider design. Lower peak temperatures are predicted if inlet pressure is raised from ambient to pi = 1.0 . Ways to do this should be implemented on slider bearing if higher load capacities are required.  51  Chapter 4 Streamline Upwind Petrov-Galerkin FE and Journal Bearings  4.1  Introduction  The streamline upwind Petrov-Galerkin finite element for thermo-hydrodynamic modeling is used for the analysis and design of journal bearings. Theoretical basis of the model applicable to plain journal bearings with requirements to account for cavitation effects was presented in chapter 2. A novel THD template for journal bearings in Cartesian coordinates that shows good convergence is developed.  4.2  Model Formulation  Figure 4.1(a) shows the classical two-dimensional journal geometry while figure 4.1(b) illustrates the model development whereby the bearing is unwrapped into regions that correspond to the bush, fluid and shaft. The oil groove presents a good reference point whereas attitude angle φ (fig 4.1(a)) depends on bearing load and speed among other factors. The computational region for pressure is mapped onto a rectangle, the x-z plane while the y axis represents the film thickness h as a function of angle θ from groove position. The minimum film thickness hmin is at an angle ϕ from the groove position as illustrated in figure 4.1a. Figure 4.2 shows the fluid region where film thickness h is a function of the x-coordinate. The cutting plane splits the oil groove in half. Periodic boundary conditions identify surfaces S3 and S5 as ‘mating’, representing the circular nature of the template.  52  The film thickness h is determined from equation 2.39. The film thickness h(θ) from equation 2.39 can be transformed to h(x) from the relation h(x)= h(θ)/(πD).  ϕ  Oil supply hole Oil groove  xb R1  yb ys W e  xs  θ  R3  R2 h  ω z x  fxs y  fys  Figure 4.1(a) Schematic representation of journal bearing.  Fluid region (exaggerated)  Shaft  Bush  4.1(b) Journal bearing unwrapped. The cross-sectional areas of the two regions are equal  53  4.1(c) Isometric view of finite element mesh of the model  Oil grooves  bush  Side surface S5 End surface S4  C2 C3  Y Z X  C4 x  C1  Mid-plane surface S1  ps bc’s  Journal (lower surface)  Side surface S3  Figure 4.2 Cross section of the film region. Note that the origin is on line C3 and surfaces S3 and S5 are defined on C4 and C3 respectively.  54  Film rupture and reformation effects have to be considered in journal bearing computations. Cavitation algorithms fall into two main categories, the Swift-Stieber (Reynolds) boundary condition which states that ∂p/∂ϕ=0 at ϕ=ϕcav, where ϕcav is the rupture interface, and the mass conserving algorithms by Elrod [28] and Kumar and Booker [29,30] among others. The Swift-Stieber condition has gained wide acceptance for steady state stationary loaded journal bearings [55] even though it does not work as well at the reformation boundary, will be used.  4.3  Simulation Results  4.3.1 Mesh Density, Convergence and Performance  The experimental data from Dowson [59] represents one of the landmark work in plain journal bearings and is often used to validate many journal bearing models. Initial model validation was performed on this data.  Input parameters and thermocouple setup  information is given in table D.1 and section D1.1. The same problem is solved analytically by Khonsari [68] using the Boyd and Raimondi tables and Yang and Khonsari’s [55] proposed empirical equations. No data on the bush dimensions and geometry of the oil supply groove thus a convective heat transfer coefficient for the bush/ambient interface of 100W/m2.K is used and found to give results close to experimental. The published peak pressure of 3.5 MPa obtained by Dowson [59] is used in convergence analysis. Table 4.1 shows the effect of varying grid size from N=10x10 to N=80x80 for both linear and quadratic elements. The sequence N=10x10, 20x20 …is used since it is recommended in numerical analysis to halve the step size in a convergence study. As expected quadratic elements converge faster although they are  55  computationally more expensive in terms of time. Figure 4.4 shows the convergence performance of the model while figure 4.5 presents the predictions for each mesh compared to the experimental value of 3.5 MPA.  Table 4.1 Convergence of Peak Pressure Predictions from Experimental Value N  Pressure Prediction  Exp. Value  Linear Elements: Devation = (Pexp-Ppred)/Ppred  Quadratic Elements: Devation = (Pexp-Ppred)/Ppred  3.50E+06  0.3128  0.2152  Linear Elements  Quadrilateral Elements  10x10  2.405E+06  2.747E+06  20x20  3.074E+06  3.194E+06  0.1217  0.0876  40x40  3.418E+06  3.453E+06  0.0235  0.0134  80x80  3.468E+06  3.476E+06  0.0093  0.0067  Percentage Divation vs No. of Elements N 35.0 Linear Quad  Percent Deviation  30.0 25.0 20.0 15.0 10.0 5.0 0.0 10  20  30  40 50 60 No. of Elements N  70  Figure 4.3 Convergence analysis on Dowson’s data  56  80  90  Pressure Predictions vs No. of Elements N 4.0 3.5  Pressure, Mpa  3.0 Pressure@N=80x80 2.5  Pressure@N=10x10  2.0  Pressure@N=20x20  1.5  Pressure@N=40x40  1.0 0.5 0.0 0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  Figure 4.4 Peak pressure prediction for Dowson’s experimental value of 3.5 MPa  Quadratic element mesh achieves a deviation of only 1.3% with a 40x40 mesh size while an 80x80 mesh is required for linear ones. Computation times as a function of mesh density for linear elements are shown in table 4.2 (using a Pentium IV 3.25 GHz 1 GB RAM). These are for a single iteration in figure 2.5. Considering that each cycle required at least three iterations for convergence and that the process is repeated at Ns=2 rev/s intervals in the optimization model, computation time can be considerable. A target deviation of 1% was used which coincides with a mesh of N=80x80 elements (0.9%).  Table 4.2 Computation Times for One Iteration  Computation Times (Sec) N 10x10 20x20 40x40 80x80  Mesh  Computing  Total  0.0313 0.0781 0.2812 1.0156  0.0156 0.0469 0.4219 1.9219  0.0469 0.125 0.7031 2.9375  57  Model performance in predicting peak temperatures for Dowson’s data is shown figure 4.5 (experimental value is 47°C). In order to compare predicted temperatures to experimental outlet drain temperatures and Khonsari’s [68] analytical computations, nodal temperatures are summed and divided by the number of nodes. Khonsari performs six iterations in his computations as shown in figure 4.5 which also indicates Dowson’s experimental value and the proposed SUPG model. The rapid convergence of the model is evident; only three iterations produce stability in predictions. The average temperature prediction of 48.4˚C from the model is 1.4˚C higher than experimental values; the peak temperature prediction of 51.2˚C compares favorably with 51.6˚C from Dowson. The peak temperature is located at the mid-plane of the journal. Table D.1 (Appendix) has input parameters for the bearing while table 4.3 shows some of the performance parameter predictions (input parameters not available from Dowson’s data were estimated). These results are preliminary; more detailed analysis follow. SUPGThermal Predictions, Analytical Method, and Experimental Temperature 65.0 T_measured T_model T_analyt  60.0 55.0  Temperature (C)  50.0 45.0  Exp. temperature  40.0 35.0  Oil inlet  30.0 25.0 20.0 1  2  3  4  5  6  7  No. of Iterations  Figure 4.5 Model, experimental, and analytical predictions, Dowson’s data.  58  Table 4.3 Model predictions, analytical [68] and Dowson’s [59] experimental data. Analytical  SUPG  Exp.  51.0 3.60 24.2 67.4 3.42  51.0 3.46 26.1 65.7 3.41  57.0 3.50 26.2 64.4 3.39  Attitude angle ϕ° Peak Pressure (Mpa) Side Leakage cm3/s Friction Force, N Shaft Torque, N.m  4.3.2  Pressure and Pressure Gradient Predictions  The experimental data from Boncompain et al. [15] and Costa et al. [69] are used to perform further validation of the model. The top surface is stationary (bush) and the bottom is the shaft in figure 4.1(b); figure 4.1(c) shows the mesh model. Input parameters for the bearings appear in Table D.2 and D.3 respectively. The effect of pressure gradients ∂p/∂x and ∂p/∂z on the energy equation 2.12 is through the compressibility term (second term on right hand side). This term is retained here even though most researchers ignore it due to its relatively small contribution [1,2,55,74]. The behavior of the two gradients near and in the rupture zone form a critical part of the discussion to follow, hence the need to retain them. Table 4.4 shows the effect of including the compressibility terms on temperature predictions. Table 4.4 Effect of pressure gradients ∂p/∂x and ∂p/∂y on temperature rise ΔT Data  Ti, °C  Dowson Boncompain Costa  36.8 40 35  ΔT, With °C  ΔT Without, °C  11.39 11.77 19.07  11.48 12.13 19.19  59  Difference, °C 0.09 0.36 0.13  Model convergence performed on pressure gradient dp/dz for the experimental data from Boncompain et al. [15] is presented in figure 4.6. The bearing analyzed has aspect ratio L/B=0.8 which is short and pressure gradient dp/dz falls rapidly from the bearing mid0.5 are due to the increased effect of film  plane. Lower values of the gradients at  pressure in the narrowing (converging) part of the bearing. Pressure predictions from the proposed model and the rupture/reformation interfaces are shown in figure 4.7 for Boncompain’s [15] bearing. The rupture zone is in the range 227 0 < θ < 332 0 which is in good agreement with Boncompain’s results (figure 4.7b). Figure 4.8a presents pressure gradient ∂p/∂x as a function of Z at several locations around the bearing. Boxed numbers indicate position along the X coordinate. Figure 4.8b compares ∂p/∂x and ∂p/∂z along the x axis and reveals that the magnitude of ∂p/∂x is much greater than ∂p/∂z, which is zero except for the rupture zone which consists of finger-like fluid/gas regions. It appears ∂p/∂x has a contribution on the compressibility term of the energy equation since it is significantly larger than ∂p/∂z although this contribution is negligible as indicated in Table 4.3.  4.3.3 Pressure Gradients and Effect on Dissipation  Of more significance is the position of the rupture interface and how it affects the gradients. The boxed numbers on the ordinate of figures 4.8a and 4.8b show the position along the x axis for Boncompain’s [15] input data (Table D2 in appendix). The magnitude of ∂p/∂x increases along the x direction, peaks at X = 0.5 , changes sign with a minimum at X = 0.592 which is located just before the rupture interface and becomes zero at the reformation interface. Noting that the gradients ∂p/∂x and ∂p/∂y are multiplied  60  by the density-temperature coefficient in equation 2.12 and w is several orders of magnitude less than u [2], the contribution of these gradients to dissipation is minimal (positive values increase dissipation). Thus investigation of the dissipation phenomenon from a mathematical point of view gives an insight as to why experimental temperatures peak just after the rupture interface and decrease thereafter.  Model Convergence for gradient dp/dz vs No of Elements N 1 0 0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  dp/dz, MN/m3  ‐1 N=10  ‐2  N=20 ‐3  N=30 N=50  ‐4 ‐5 ‐6 ‐7  Figure 4.6. Convergence analysis on pressure gradient dp/dz.  61  0.9  1.0  Model Pressure Prediction, MPa 5 4.5  THD pressure  Pressure, MPa  4 3.5 3 2.5 2 1.5 1 0.5 0 0  40  80  120 160 200 240 280 320 360  Angle from oil hole, degrees  (a)  (b)  Cavitation zone  Rupture interface  Reformation interface  (c) Figure 4.7(a) Model pressure predictions at bearing mid-plane, model vs (b) Boncompain’s result [15]. (c) Pressure map showing the rupture and reformation interfaces.  62  Pressure gradient dp/dx  150 4 100 3 5 50 6  1, 2,10,11  0  dp/dx, MN/m3  9  0.00  0.10  0.20  0.30  0.40  0.50  7 -50 12 -100  1. dpdx@X=0.11 2. dpdx@X=0.21 3. dpdx@X=0.46 4. dpdx@X=0.5 5. dpdx@X=0.53 6. dpdx@X=0.54 7. dpdx@X=0.55 8. dpdx@X=0.61 9. dpdx@X=0.65 10. dpdx@X=0.92 11. dpdx@X=0.99 12. dpdx@X=0.63  -150  -200 8 -250  (a)  Gradients dp/dx and dp/dz, MN/m3  Pressure Gradient Variations 150 100  dp/dx  dp/dz  50 0 ‐50 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  ‐100 ‐150 ‐200 ‐250  (b) Figure 4.8(a) Pressure gradient ∂p/∂x variation along the bearing x coordinate, Costa et al. [69]. (b) Pressure gradients ∂p/∂x, ∂p/∂z variation along the x axis.  4.3.4 Dissipation Source Terms  The dissipation source terms (last two terms in equation 2.12) in the energy equation are evaluated for their contribution to the thermal profile. Figure 4.9 shows the values of the  63  gradient ∂u/∂y (normalized by dividing by the peak value) as a function of x . Peak values occur at the rupture interface 1 which is just before peak temperatures. This is an important observation as it correctly indicates the region of maximum viscous dissipation without specifying thermal boundary conditions at this interface as done by Boncompain [15] and other researchers [6]. Velocity gradient ∂w/∂y is shown in figure 4.10 (normalized by dividing by peak ∂u/∂y and multiplying by 100). An order of magnitude analysis reveals this gradient varies between -0.38 percent to +0.18 percent of ∂u/∂y making its contribution small and hence justifies ignoring the term [2]. The peak value occurs at X = 0.98 which is close to the oil groove. Discontinuities (figure 4.10) appear to be caused by the groove/cavity interface and mixing zone phenomenon.  64  (a)  (b) Figure 4.9(a) Velocity gradient ∂u/∂y as a function of X . (b) Magnified view at the peak gradients point.  65  Figure 4.10 Velocity gradient ∂w/∂y as a function of X .  4.3.5 Thermal Profiles  One benefit of the model is the detailed thermal profile of the film. Figure 4.11(a) shows the lower surface is in contact with the shaft and the top with the bush. The bottom surface is assumed to be at a constant temperature whereas the bush contact surface varies significantly. As mentioned above temperature measurements are normally taken on the bush or shaft; knowledge of the fluid profile is sketchy at best. The isotherms in 66  the axial (z) direction are mainly horizontal indicating no thermal gradients confirming what other researchers noted [2,67]. Their density depends on the number of transverse divisions considered (ten in the model). The variation in temperature across the film thickness is apparent due to the mixing of the fresh oil with circulating oil. Comparison of proposed model predictions (figure 4.11(a)) and figure 4.11(b) (figure 7, adopted from Boncompain et al. [15]) at the minimum film thickness position (indicated by the 0 centerline) and position x ≈ 0.1(36 ) in figure 4.11(a) shows similarities. The circular  peak just after hmin in the figure is also very similar to the circular peak just after hmin in figure 4.11(b). A peak temperature of 53.5°C occurs at an angle of 232° from the oil supply groove (position (0.65,0,0.12)). The geometry of the oil supply hole and cavity are not given in [15] which makes setting boundary conditions difficult. Conservation of mass flow at the inlet was used to determine oil temperature after mixing. Experimental data from Costa et al. [69] and Kucinschi et al. [27] are analyzed next. Input parameters for the bearings are shown in Table D.3. Temperature profiles from the proposed model are presented in Figure 4.12(a) and (b), including the inlet hole and cavity as a function of inlet boundary conditions (note that the x-axis compressed due to large aspect ratio between this axis and the z-axis). The proposed model predicts Costa’s data well. Note that Costa et al. [69] present data and no model in their paper. In figure 4.12(a) the supply hole and cavity are assumed to be at Ti whereas in 4.12(b) it is assumed to be at Ti. The latter assumption results in higher cavity temperatures. A peak temperature of 62.7°C at position x ≈ 0.70 (252°) is in agreement with 62°C (figure 5 in Costa [69]) although the position is slightly lagging. The cavity thermal profile depends on boundary conditions imposed on the supply hole and cavity. In figure 4.12(a) both the  67  supply hole and cavity are assumed to be at temperature Ti whereas in figure 4.12(b) the entrance to the supply hole is at Ti. Clearly a mixing phenomenon is evident in the latter case which appears closer to reality. The type of boundary condition used plays an important role in model predictions. Kosasih and Tieu [5] assumed that most parts of supply region are filled with re-circulating oil and most of the injected cold oil does not immediately enter the oil film, which tends to support the case in figure 4.10(b). The influence of the inlet oil at a lower temperature Ti on the cavity oil dynamics is also clearly evident. Figure 4.13(a) shows the bush upper thermal profile while in figure 4.13(b) the bush lower surface is exposed to show how the temperature varies along the axial direction (note that both coordinates y, z are reported in real figures and not in non-dimensional form). Coordinate z is the actual bearing half width dimension and y is magnified 4000 times since it is very small compared to the other two). Model results of the bush temperatures are in agreement with Costa’s et al. [69] experimental results. The peak temperature of 62.7°C occurs at position x ≈ 0.70 . As expected the outer edge ( z = 0.04 ) has lower surface temperatures than the mid-plane due to convection to the ambient and leakage flow (the top isotherm of 51.3°C does not get to the outer surface).  4.5.3.1 Standard Error and Goodness of Fit  To validate the model further, Costa’s experimental temperature data at the mid-plane of the bush are compared to model predictions in figure 4.14 and Table 4.3. Standard error  σe is used to estimate the error ∑  ∑  (4.3)  68  where n is the number of points in each series, i the point number in series s, yis the data value of series s and the ith point, and ny the total number of values in all series. The standard error of a method of measurement or estimation is the estimated standard deviation of the error in that method. Specifically, it estimates the standard deviation of the difference between the measured or estimated values and the true values. Notice that the true value of the standard deviation is usually unknown and the use of the term standard error carries with it the idea that an estimate of this unknown quantity is being  used. Looking at figure 4.14 model predictions at θ=-30° has the most points outside the standard error estimate of both the data and predictions (9 points) whereas all points at  θ=0° are within and 2 points are outside at θ=+30°. Table 4.5 summarizes predictions from the three positions (θ=-30°,0,+30) including the square of the Pearson correlation coefficient r2. A possible explanation for the high error for θ=-30° is that this groove location might be in the cavitation zone and so applied boundary conditions are not applicable. Table 4.5 Standard Error and r2 correlation for Temperature Predictions θ=-30 θ=0 θ=+30 Standard Error r2  3.18  2.12  2.66  0.796  0.939  0.921  Note: r2 represents the square of the Pearson product moment correlation coefficient through data points. As explained in section 2.4.2 the inlet film temperature Ti can be determined using 1  conservation of energy in the groove  where Ts, Tr are the  supply and re-circulating oil temperatures respectively and Qr, Qin the re-circulating and oil inflow. Maximum model prediction deviation from experimental values is 0.09 which  69  occurs in the region after oil mixing prior to the minimum film thickness at groove location θ=-30°. A probable cause of the higher error here is that the groove is probably in the rupture zone. This indicates a need for a more accurate modeling of the mixing phenomenon for this case. Peak temperature predictions are within 0.03 in each case. It has to be emphasized here that unlike other models in the literature [20,27], the model as presented is dynamic in that groove location can be changed and results may be easily obtained. Appendix E explores predictions from other models and shows the superior performance of SUPG in steady state THD lubrication. Peak experimental temperatures as a function of supply pressure for a 2 kN load and 50 rev/s speed are shown in figure 4.15. Boxed values indicate experimental temperatures while standard error bars are also included. Deviation of model predictions is small as shown in Table 4.7. Experimental values differ from model prediction by a maximum of 1.6 percent.  70  (a)  (b) Figure 4.11 (a) Thermal profile from Boncompain et al. [15] data. (b) Figure 7 from Boncompain. The geometry of the oil supply hole and cavity are not given which makes setting boundary conditions difficult.  71  1 08  1 0.8 0.6  (b) Figure 4.12(a and b) Thermal profile in the inlet, cavity and bearing, Costa.  72  Y+  0.4 0.2 0  Y+  06 04 02 0  (a)  (a)  (b) Figure 4.13 Bush thermal profile, (a) top and (b) bottom view  73  Model Prediction vs Exp. temperatures, θ=-30 65  Temperature, C  60 55 50 45 SUPG Exp. data  40 35 30 0  30  60  90  120  150  180  210  240  270  300  330  360  330  360  Angle, degrees  (a) Model Prediction vs Exp. temperatures, θ=0 70 65  Temperature, C  60 55 50 SUPG Exp. data  45 40 35 30 0  30  60  90  120  150 180 210 240 Angle, degrees  270  300  (b) Model Prediction vs Exp. temperatures, θ=+30 65  Temperature, C  60 55 50 45 40  Exp. data SUPG  35 30 0  30  60  90  120 150 180 210 240 270 300 330 360 Angle, degrees  (c) Figure 4.14 Thermal predictions at the bush/oil interface – Costa’s data. Input parameters given in Table D3.  74  Table 4.6 Predicted versus experimental temperatures, Costa’s data. Exp. Values [69]  x 0.0000 0.0080 0.0109 0.0393 0.0654 0.0916 0.1178 0.1309 0.1440 0.1571 0.1702 0.1833 0.1963 0.2182 0.2313 0.2443 0.2662 0.3033 0.3107 0.3142  Model Predictions θ =-30  Angle (Degrees)  T(-30) Fillon  T(0) Fillon  0 9 13 45 75 105 135 150 165 180 195 210 225 250 265 280 305 348 356 360  35.0 35.0 35.5 39.0 41.8 42.0 43.0 44.0 45.0 46.5 47.0 49.0 53.0 55.0 58.0 59.0 58.0 46.0 35.0 35.0  35.0 35.0 35.5 39.0 41.0 42.0 44.0 46.0 47.5 49.0 51.0 54.0 58.0 60.0 61.0 62.0 58.0 46.0 35.0 35.0  θ =0  Difference %  θ =+30  T(+30) Fillon  Diff. % Diff. % Diff. % (-30) (0) (+30)  35.0 35.0 35.0 35.0 35.0 35.0 35.5 36.4 36.3 37.0 42.6 42.9 38.0 46.0 45.5 39.0 46.7 46.7 44.0 47.8 48.9 49.0 48.9 51.1 50.0 50.0 52.8 52.0 51.7 54.4 54.0 52.3 55.3 58.0 54.5 56.2 60.0 57.1 57.2 60.5 57.8 57.9 60.0 57.7 57.7 58.0 56.9 56.9 56.0 53.2 53.2 40.0 42.0 41.8 35.0 35.0 35.0 35.0 35.0 35.0 Input Parameters Bearing load, kN  35.0 35.0 36.3 41.1 42.2 43.3 48.8 53.9 54.6 55.4 56.2 57.1 58.1 58.7 58.6 57.8 53.9 39.8 35.0 35.0  Bearing speed, Ns rev/s Supply pressure, Ps, kPa Inlet temperature, C Inlet viscosity, Pa.s All other input parameters same as Table D3  75  0.0 0.0 2.5 8.5 9.1 10.0 10.0 10.0 10.0 10.0 10.0 10.0 7.2 4.8 0.5 3.7 9.0 9.5 0.0 0.0  0.0 0.0 2.2 9.1 9.9 10.0 10.0 10.0 10.0 9.9 7.8 3.9 1.4 3.6 5.7 9.0 9.0 10.0 0.0 0.0  0.0 0.0 2.2 10.0 10.0 9.9 9.8 9.1 8.4 6.1 3.9 1.6 3.3 3.1 2.4 0.3 3.9 0.5 0.0 0.0  8.0 50.0 300.0 35.0 0.0  Temperature as function of supply pressure 65 64.5 64  Peak Temperature, C  63.5  64.4  63.4  63 T_model  62.5  Data  62  63.9  61.5 63.2  61 62.5  60.5 60 59.5 0  50  100  150 200 Supply Pressure, kPa  250  300  350  Figure 4.15 Predicted versus experimental peak temperatures as function of oil supply pressure, Costa et al. [69]. The y-axis starts at 60.5°C hence the apparent large difference between prediction and experimental data.  Table 4.7 Predicted versus experimental temperatures as function of supply pressure. Model Exp. Predictions, Values, °C (Angle) °C 63.5 (300°) 64.4 63.0 63.2 (298°) 63.2 62 (293°) 62.0 61.4 (288°) 61.2 60.7 (279°) Load, kN Groove location θ Bearing Speed, Ns rev/s Oil supply temperature, C  Supply Pressure, kPa 50 70 140 200 300  76  Difference % 1.40 0.32 1.90 0.97 0.82 2 (+30) 50 35  4.3.6 Pressure Simulations of Costa’s Data  Pressure predictions presented some challenges, the main one being how to take the effect of supply pressured into account. It appears from Costa et al. [69] that supply pressure is insignificant at groove locations (theta=-30 and 0 degrees). However, an initial pressure of 200 kPa is detected at groove location theta=+30. Figure 4.16 presents model predictions versus experimental pressures. Standard error bars show predictions to be within the ‘confidence’ range.  77  Pressure Predictions vs Exp. data, θ=-30 1800 1600  P_Model  Pressure, kPa  1400  Exp. data  1200 1000 800 600 400 200 0 0  30  60  90  120  150  180  210  240  270  300  330  360  Angle, degrees  (a) Pressure Model Prediction Vs Measured, Theta=0 2000 1800  Pressure, kPa  1600 P_model P_measured  1400 1200 1000 800 600 400 200 0 0  30  60  90  120  150  180  210  240  270  300  330  360  Angle, Degrees  (b) Pressure Prediction vs Model, Theta=+30 2000 1800  Pressure, kPa  1600 1400  P_exp P model  1200 1000 800 600 400 200 0 0  30  60  90 120 150 180 210 240 270 300 330 360  Angle, degrees  (c) Figure 4.16 Pressure prediction as a function of groove location, Costa et al. [69]. Applied load: 6 kN, speed: 50 rev/s.  78  4.3.7 Flow Predictions  Flow predictions are compared to experimental data [69] at a bearing load of 8 kN, 50 rev/s speed, and an inlet angle of -30˚ in figure 4.17. An analytical prediction from Martin’s formula [70] (see Appendix B) is also included for comparison. The model predictions are within 2.8 percent of measured values while Martin’s formula gives an accuracy of 11.2 percent. Standard error bars are included for both models. More detailed analysis of flow predictions is presented in chapter 5.  Flow Predictions at W=8kN, Theta=-30 2.8 2.6 2.4  Flow, L/min  2.2 2 1.8 Martin,W=8kN,Theta=-30 Measured Flow_model  1.6 1.4 1.2 1 0  30  60  90  120  150  180  210  240  270  300  330  Supply Pressure, kPa  Figure 4.17 Flow predictions: model, experimental data and Martin’s formula [70].  4.4  Chapter Summary and Conclusions  An effective and systematic procedure to model the thermo-hydrodynamic phenomenon is presented. The cylindrical coordinates are transformed to Cartesian ones with the resultant ‘bearing template’ used to analyze plain journal bearings. This transformation  79  makes bearing analysis easier as shown by the rapid convergence of the model. On average only two iterations are required to produce stability in thermal predictions. This difference between model and experimental data is consistently less than 10 percent. No special treatment of the rupture interface is required in this model due to the high stability of the streamline upwind Petrov-Galerkin method. The velocity gradient ∂u/∂y is maximum at the rupture interface while ∂w/∂y, which is significantly less, has its lowest values after the reformation interface (no side flow in this region). The model also shows the high influence boundary conditions at the oil inlet/groove have on the results. More investigations need to be carried out to determine the most realistic boundary conditions in the mixing zone.  Improvement on the mixing model by  implementing conservation of energy in the groove led to improved predictions shown in Appendix E. Pressure simulations are more difficult due to a lack of knowledge about the effect of supply pressure on inlet boundary conditions. Peak pressure predictions are however, reasonably good. Flow predictions from the model show a maximum deviation of 2.8 per cent as compared to predictions from Martin’s formula. The model is thus a significant improvement on methods currently available. Further verification of SUPG appears in Appendix E where the method is shown to be effective in convection-diffusion problems such as the one analyzed by Kelly et al. [71]. Improvements to the inlet boundary condition modeling by implementing groove mixing using conservation of energy resulted in improved temperature predictions for Mitsui’s [72] and Pierre’s et al.[73] data (Appendix E). Prediction accuracy of the SUPG model  80  can be improved significantly with better groove mixing computations which are complex due to groove size and geometry.  81  Chapter 5 Formulation of Empirically Derived Model and Simulation Results  5.1 Introduction  The SUPG finite element model developed in the previous chapters is used to validate optimization equations that are introduced here to complete the integrated design process for journal bearings. Empirical equations for temperature increase ΔT, leakage flow QL and power loss PL are formulated from first principles. Coefficients for these equations are determined by simulating measured bearing data in the literature. For a given set of design variables the FE model is run at speeds from zero to 200 rev/s and the three parameters are determined at each speed. Figure 5.1 gives details of how the optimization model is structured. The main focus in this work is to use the developed THD model to validate a proposed optimization scheme for medium range (length to diameter ratio 0.4<λ<1.5) medium speed (0<Ns<200 rev/s) plain journal bearings. Journal bearing design and optimization involves choosing a set of design variables that include but are not limited to radial clearance C, length to diameter ratio λ, viscosity μ, oil supply pressure ps for optimum performance and groove location θ. The goal is to find a set of design variables that minimizes the objective function. The objective function may include temperature increase ΔT, side leakage QL, and power loss PL. Derivation of equations for temperature rise ΔT, side leakage QL, power loss PL, and constraint equations for maximum pressure pmax, whirl onset speed ωcr, and oil supply pressure ps will be considered. It should be  emphasized that selecting design variables may be influenced by additional factors such  82  as journal strength/stiffness requirements for minimum misalignment. Thus the minimum journal size could also be determined by consideration of shaft deflection and strength which is not considered here. A formal definition of the optimization model and its formulation is presented in chapter 5. Multi-objective function equations with linear and non-linear constraints form the basis of the optimization models developed. Two solution techniques, namely the method of assigning weighting/scaling factors and using Pareto optimal fronts will be explored. The targeted range of aspect (L/D) ratios is 0.4<λ<1.5 and a speed range (0<Ns<200 rev/s). Although Hirani and Suh [54] propose an optimization model for this range of bearings, their model requires complex numerical integration for power loss PL and use of mass conserving finite difference equations to solve for angle θ  around the bearing  circumference. Their model also involves reference to tables (table 1 in reference 38) which makes their model even more cumbersome to implement as a flow prediction program. Optimization equations are derived from first principles and empirical validation of coefficients is achieved using the bearing template developed in the previous chapter as well as proving these equations using data from literature. Convergence and simulation times for both meshing and computations was addressed in section 4.3. As noted each iteration took about 2.9 seconds and at least three iterations were required for convergence for a total of 8.7 seconds. For optimization computations which required obtaining operating characteristics at Ns=2 rev/s intervals, a typical cycle time for Costa’s data (Tbale D3 in appendices) is analyzed in table 5.1.  83  Table 5.1 Computation Times for Optimization Model Time Stepping  Task Meshing Time for convergence at each Ns Total time for 0<Ns<200 rev/s Time required for SQP algorithm Time for seven attempts with SQP Time required for hybrid algorithm  Time (Seconds) 4.1 26.6 2510.5 5.0 35.0 15.0  The total time of 2510 seconds translates to approximately 42 minutes for the iterations. The computer specifications are Pentium IV clock speed of 3.25 GHz and 1 GB RAM. The genetic algorithm generally required more time than SQP, however the problem with SQP is the result depends on the initial guess. If the initial guess is in a region with a local minimum the algorithm gives erroneous results. On average seven attempts were made to get reasonable results with SQP. The SUPG finite element method predicts temperature increase from given input parameters and further validates the coefficients in equation 5.7. Table 5.2 presents the coefficients.  84  Start with initial speed Ns  Run SUPG FE model to compute ΔT1, QL1, and PL1  Predict ΔT2, QL2, and PL2 using empirically derived model Increase speed Ns at 2 rev/s intervals up to 200 rev/s  Readjust coefficients for empirical model  No  Are predictions (ΔT2, QL2, PL2 ) within specified tolerances to (ΔT1, QL1, PL1) and data in literature?  Yes Empirical model validated at current speed Ns  Determine optimized design variables using Hybrid Technique  End  Figure 5.1 Flowchart for Optimization Model  85  5.2  Governing Equations  5.2.1 Temperature Increase  The journal bearing shown in figure 4.1(a) is used for reference. The method proposed follows Jang and Khonsari’s approach [55] whereby an initial estimate of the Sommerfeld number S is made  S=  μ eff N s DL ⎛ R ⎞ 2 W  ⎜ ⎟ ⎝C ⎠  (5.1)  where the oil viscosity at inlet conditions is taken as the effective viscosity μeff. Eccentricity ratio ε is determined to be a function of Sommerfeld number and length to diameter ratio λ, ε=ε(S,λ) using the approximation  ε = −0.36667 +  0.49737  λ  − 0.20986 ln (S )  0.5 ≤ λ ≤ 1.5  (5.2)  A closed form of the energy equation can be derived by assuming adiabatic conditions at the oil/bush interface ( ∂T  ∫  2π  0  ∂T dx = 0 ∂y  ∂y  = 0 at y = 0 ) and  and  T = Tshaft  at y = 1  (5.3)  while isothermal conditions apply at the shaft/oil interface ( T = Tshaft at y = 1 ), the socalled ISOADI (isothermal shaft, adiabatic interface) condition resulting in the dimensionless energy equation [74], as presented in equation 2.14b. Two key parameters  κ1 and κ2 that appear in equation 2.14b are referred to as temperature rise parameters. The first temperature rise parameter κ1 is associated with viscous dissipation while κ2 relates oil properties and velocity. The energy equation is solved for temperature increase ΔT(X); it is evident ΔT(X) is a function of eccentricity ratio, κ1 and κ2  86  ΔT(X) = f(ε,κ1,κ2)  (5.4)  The dimensionless temperature equation proposed by Jang [55] forms the basis for temperature increase computations ⎡ b c ⎤ Tmax = exp ⎢ a + 0.5 + 0.5 ⎥ κ1 κ2 ⎦ ⎣  (5.5)  where constants a, b, and c are empirically determined from simulating bearing temperature data in the literature [59,71,75].  Table 5.2 Values a, b, and c computation of maximum temperature Tmax. Note that the valid range of temperature rise parameters is: 0.001≤κ1≤0.5 and 0.01≤κ2≤5  Tmax = exp[a + b / κ 10.5 + c / κ 20.5 ] Eccentricity ratio 0.1<ε<0.3 0.3≤ε≤0.7 0.7<ε≤0.9  Constants a 1.58105 1.5917 1.63973  b -0.14250 -0.13962 -0.10627  c -1.04050 -1.24999 -1.45216  These constants are functions of the eccentricity ratio. Temperature predictions with the constants in table 5.2 are compared to the finite element and other models in section 5.4 and shown to give accurate results. The maximum temperature can be determined from the relation Tmax = β (T − Ti ) once Tmax is known. Maximum shaft temperature is determined in a similar manner  Tshaft  ⎡ b c ⎤ = ⎢a + 0.5 + 0.5 ⎥ κ1 κ2 ⎦ ⎣  −1  (5.6)  Effective viscosity μeff is updated using the shaft temperature Tshaft and equation 2.26. The Sommerfeld number and eccentricity ratio are updated using the new viscosity  87  values; maximum temperature and shaft temperature are re-evaluated using the updated values. The iteration process is terminated when convergence occurs to a predetermined value.  5.2.2  Side Leakage  A key difference with other models is in the empirical relation for the leakage flow QL . Khonsari’s model [55] utilizes QL = QL =  π 4  π 2  N s CDL QL  while Hashimoto [34] uses  N s CD 2 ε for short bearings. Consideration of where the hydrodynamic film starts  along the bearing circumference is important [68], however the curve fitted empirical model developed here appears to give good results. The non-dimensional flow QL is expressed as a function of ε and λ as proposed by Martin [70] and used by Jang and Khonsari [55]  QL = exp[0.20 − 0.500λ1.5 + 1.05 ln(ε )]  0.5≤λ≤1.5; 0.1≤ε≤0.9  (5.7)  where the coefficients were determined empirically by Gadala and Zengeya [76] who simulated pressure-fed bearings in the literature and performing curve-fitting to the data given that QL =f(ε,λ). Model predictions are compared to those from Khonsari’s model in simulating Costa’s data in section 4.4. Predictions from Khonsari’s model are not as accurate since they don’t specifically consider supply pressure ps. The proposed model was also validated using measured leakage data from literature [59,69,75] and compared with predictions from Song [36], Jang and Khonsari [55], and the finite element model.  88  5.2.3 Power Loss Determination  Power loss forms an important consideration in journal bearing design. Frictional losses in machines result in wasted energy and the generation of heat which affects the life of materials including the lubricant. A significant portion of total fluid film bearing losses in high-speed turbomachinery may be consumed simply in overcoming frictional losses. At surface speeds below approximately 51 m/s, these losses are generally less than 10 percent of total power input to the bearing [77]. At higher velocities, however, they rapidly with surface speed and can reach 25 to 50 percent of the total bearing loss in large turbine-generators and related machinery. Power loss is determined from ⎛ μ eff U h ∂P ⎞ ⎟URdθdz PL = ∫∫ ⎜⎜ + 2 R ⋅ ∂θ ⎟⎠ ⎝ h  (5.8)  where the effective viscosity is determined at the current iteration. Hirani and Suh [38] use a mass conserving finite difference scheme to solve for the power loss. The drawback in this approach is that it involves much more computational effort than proposed. The approach in emanates from the established fact the term (R C) f is a function of λ and ε. Power loss is determined from  ⎛ ⎞ R 2.73622 − 4.97838 ε ⎟⎟ f = exp⎜⎜1.9994 + C λ ⎝ ⎠  (5.9a)  Ffr=fW  (5. 9b)  PL=Ffr(2πRNs)  (5.9c)  where equation 5.9a is determined empirically following Jang and Khonsari’s approach [55]. The frictional force Ffr is determined once f is known, and W the load on the bearing.  89  5.2.4  Pressure Computations  Pressure prediction for the optimization equations are implemented by integrating the Reynolds equation to get a non-dimensional form [55] p − ps P= μi N s  ⎛R⎞ ⎜ ⎟ ⎝C ⎠  2  (5.10)  with the empirical formula proposed by Jang and Khonsari [55] used to evaluate nondimensional pressure P  0.39491 ⎡ ⎤ + 4.77721ε 3 ⎥ Pmax = exp⎢2.25062 − λ ⎣ ⎦  (5.11)  The onset of whirl is important due to the instabilities it creates in journal bearings. However, for speeds considered here (Ns<200 rev/s) is was determined not to be a major problem even for the upper range. Hashimoto [34] gives a brief account on the effects of whirl especially for high-speed, short bearings. Whirl velocity is determined from the Routh-Hurwitz stability criterion and found to be a function of eccentricity ratio ε. An empirical expression is presented  ωcr =  g {0.0584 exp(6.99ε 2.07 ) − 1.318ε + 2.873} C  (5.12)  where ωcr is the critical whirl speed and g the gravitation acceleration. 5.2.5  Formulation of the Optimization Problem  The first part of the work investigated minimization of temperature increase ΔT(X) and side leakage QL by combining the two objective functions F1 and F2 into a multiobjective function F(X) using weighting and scaling parameters. This works well if the relative importance of each functional is known apriori. Hashimoto [34,35] uses this method in analyzing high-speed, short journal bearings. The concept of axiomatic design  90  [39] helps determine the independence of objective functions, and as explained by Hirani [40,41] ΔT (F1) and QL (F2) are not independent of each other as ΔT is a function of W and QL which means F1=func(W, F2). However objective functions W (F1) and QL (F2) are making them the preferred choice. In addition to the method of weighting functions Pareto optimal fronts at several common speeds are produced as an alternative optimization technique that requires no weighting or scaling parameters. Pareto optimal fronts are sets of equally efficient or non-inferior alternative design vectors. The design vector X is considered non-inferior if the following conditions are fulfilled [41] Wi≤Wj and Qi<Qj or  (5.13a)  Wi<Wj and Qi≤Qj  (5.13b)  Four design variables are chosen; namely the radial clearance ratio C, length to diameter ratio λ, oil viscosity μ and supply pressure ps. The design variable set can be written in a vector X, X T = (C , λ , μ , p s ) . Supply pressure is specifically included as emphasis is on pressure-fed bearings and variation of this variable will be closely monitored in the results. State variables vary with the given operating conditions of the bearing such as the load W, rotational speed Ns and ambient temperature To. These include the eccentricity ratio ε,  film pressure p, film temperature T (K), side leakage QL (m3/s), and whirl onset speed  ωcr. The given choice of design variables affects the state variables. It should be emphasized that selection of design variables may be influenced by additional factors such as journal strength/stiffness requirements for minimum misalignment. Thus the minimum journal size could also be determined by shaft deflection which is not considered here.  91  The multi-objective optimization problem can be formally stated as [42] Find X to minimize F(X) subject to the constraints gi(X) ≤ 0 where F(X) = α1β1PL + α2β2QL  (5.14a)  and gi(X) ≤ 0,  (5.14b)  with  g1 = ha – C{1 - ε(X)},  g2 = ΔT(X) – ΔTa,  g3 = ω - ωcr,  g4 = pmax – pa  non-linear constraints and Cmin ≤ C ≤ Cmax,  λmin ≤ λ ≤ λmax, μmin ≤ μ ≤ μmax, ps, min ≤ ps ≤ ps, max  (5.14c)  are lower and upper bounds linear constraints for C, λ, μ and ps respectively. Constants  α1 and α2 are scaling while β1 and β2 fulfill a weighting role in equation 5.14a. Note that these are not necessary if the problem is solved using Pareto optimal fronts.  5.2.6  Implementation in Matlab  The program structure follows the steps given in Appendix C and is implemented in Matlab [51,78] with the input data in table D.5. The hybrid scheme invokes the genetic algorithm (GA) for rapid convergence to the global minimum region and SQP once this region is identified for accurate determination of the minimum point. Equation 5.14(a) is the fitness function and design variable set X T = (C , λ , μ , p s ) comprises the set of genes. Highlights of the GA include a population of 100 individuals and a double vector population type, an initial penalty of 10 and a penalty factor of 100 for the non-linear constraints. The penalty factor increases the penalty parameter when the problem is not solved to required accuracy and constraints are not satisfied. A rank-based fitness scaling is used. Rank scales the raw scores based on the rank of each individual, rather than its  92  score. The rank of an individual is its position in the sorted scores. The rank of the fittest individual is 1, the next fittest is 2 and so on. Rank fitness scaling removes the effect of the spread of the raw scores. The selection function used is stochastic. Stochastic uniform lays out a line in which each parent corresponds to a section of the line of length proportional to its expectation. The algorithm moves along the line in steps of equal size, one step for each parent. At each step, the algorithm allocates a parent from the section it lands on. The first step is a uniform random number less than the step size. However as Hirani and Suh [38] point out it is important to maintain diversity in non-inferior solutions by including a degree of sharing. Reproduction options determine how the next generation is created. Elite individuals are guaranteed survival to the next generation; the elite count was set to 2. Crossover fraction, which specifies the fraction of the next generation other than elite individuals that are produced by crossover, was set at 0.8. The remaining individuals, other than elite individuals, in the next generation are produced by mutation. Mutation functions make small random changes in individuals in the population, which provide genetic diversity and enable the GA to search a broader space. The mutation method used is Gaussian which adds a random number to each vector entry of an individual. This random number is taken from a Gaussian distribution centered on zero. The variance of this distribution can be controlled with two parameters. The Scale parameter determines the variance at the first generation. The Shrink parameter controls how variance shrinks as generations go by. If the Shrink parameter is 0, the variance is constant. If the Shrink parameter is 1,  93  the variance shrinks to 0 linearly as the last generation is reached. Both scale and shrink parameters are set at 1. Crossover combines two individuals, or parents, to form a new individual, or child, for the next generation. The scattered option, which creates a random binary vector and then selects the genes where the vector is a 1 from the first parent, and the genes where the vector is a 0 from the second parent, and combines the genes to form the child, is chosen. Migration is the movement of individuals between subpopulations if population size is greater than one. Every so often, the best individuals from one subpopulation replace the worst individuals in another subpopulation. Migration was controlled by using the forward option where migration takes place toward the last subpopulation, that is the nth subpopulation migrates into the (n+1)th subpopulation. A migration fraction of 0.2 is used, and this controls how individuals move between subpopulations. A migration interval of 20 is set, and this controls number of generations before migration takes place. The Matlab code is included in Appendix C.  5.3  Validation of Proposed Equations  5.3.1  Effect of Supply Pressure ps  The effect of supply pressure ps on leakage flow with model predictions on Costa’s data [69] is shown in figure 5.2. Jang and Khonsari’s [55] predictions are also included for comparisons. Although Khonsari does not specifically include supply pressure in the simple expressions of Table 4, reference 55, predicted values in figure 5.2 for this model are not constant. This may be due to an indirect effect of supply pressure on eccentricity ratio ε and the way coefficients of QL are determined in the (ref. 55). Standard error and the square of the Pearson product moment correlation coefficient appear in table 5.3. It is  94  interesting to note that r2 for Khonsari’s model is better (0.999) than the model (0.998) although it is evident from the graph that the ‘goodness of fit’ is poor.  Leakage Flow Data: Costa's Bearing 70E-6 QLModel, 8kN QLKhonsari, 8kN Measured, 8kN QLModel@2kN QLKhonsari@2kN Measured@2kN  60E-6 50E-6  QL , m 3/s  40E-6 30E-6 20E-6 10E-6 000E+0 0  50  100  150  200  250  300  350  Supply Pressure, kPa  Figure 5.2 Model versus Khonsari’s predictions – Costa’s data.  Table 5.3 Standard error and square of the Pearson correlation  2 kN  5.3.2  STD Error r2  Model  Khonsari  8.87E-07 0.986  2.88E-06 0.857  Temperature Increase ΔT(X)  The equations were validated by comparing predictions to measured data from literature. Three bearings with input data in Tables D.1 (Dowson et al. [59]), D.3 (Costa et al. [69]) and D.4 (Ferron et al. [75]) are used for the analysis. Temperature increase ΔT(X) as a function of speed for Dowson’s bearing is shown in figure 5.3(a) and 5.3(b), figures  95  5.4(a) and 5.4(b) for Ferron’s data (ΔT is 18 °C at 66.7 rev/s), while figures 5.5(a) and (b) show simulations for Costa’s data. The dashed lines indicate measured data points and journal speed on all the figures. The first figure in each case has all four models (FE, proposed model, Khonsari, and Song); the second shows a smaller speed range for better clarity of measured points/values. Temperature increase ΔT(X) is determined as the leakage oil temperature while model predictions represent the average of all nodal points in the film mesh. In all cases in the figures, temperature increase is assumed to be the measured temperature of leakage oil. Standard error bars indicated on 5.3b, 5.4b and5.5b show model predictions in the 95% confidence level. In all cases the proposed model predictions are more accurate as the formulations are based on both theoretical and empirical validation for the range of bearings under consideration. The equations were derived from first principles using ISOADI and adjusted empirically using Jang and Khonsari [55] as a guide. Also included is Hashimoto/Song model [34,36] for short, high speed bearing configurations. It is clear from figures 5.3(a) and 5.4(a) that predictions from Hashimoto are inaccurate for finite and long bearings under consideration. The model is, however, included for completeness in the discussion of journal bearing optimization. A different perspective to the results is shown in Figure 5.6 where temperature increase is shown as a function of diameter and radial clearance for Dowson’s data. Measured values are indicated on the graphs. An important consideration is establishing whether the range of bearing performance is laminar and below the speed at which shaft whirl starts. This is generally determined by  96  computing the Reynolds number (Re=ρCU/μ). A sampling of data points was used where the variable set XT=(C,λ,μ,ps) for some of data is shown in Table 5.4.  Table 5.4 Local Reynolds number from a sample of the data XT in (μm,-,Pa.s,kPa)  Bearing Diameter, m  Re  (60,0.4,0.03,300) (50,0.4,0.03,300) (80,0.4,0.03,300) (40,0.6,0.03,300)  0.191 0.191 0.191 0.128  2871 3402 2284 1511  It was determined from analysis that the range of Reynolds numbers determined were in the laminar range although the second value in Table 5.4 is high. Onset of whirl occurred at Ns>195 rev/s for this particular set and higher speeds should be avoided. For Dowson’s bearing,  a  Reynolds  number  of  820  is  determined  at  Ns=200  rev/s  for  XT=(40,0.75,0.03,300) which is well below both whirl onset and turbulent flow. Costa  and Ferron bearings also operate below turbulent flow.  97  Delta T Comparisons, Dowson's Bearing 800 700 600  ΔT, C  500 DT_model DT_FE DT_Khonsari DT_Song/Hashimoto  400 300 200 100 0 10  30  50  70  90  110  130  150  170  190  210  Speed, rev/s  (a) Model, FE, Khonsari and Song/Hashimoto ΔT Comparisons, Dowson's Bearing 220 200 180 160  ΔTmesured at 25rps = 14.7 C  ΔT, C  140 120 Exp. Point DT_FE DT_model DT_Khonsari DT_Song/Hashimoto  100 80 60 40 20 0 5  10  15  20  25  30 35 Speed, rev/s  40  45  (b) Figure 5.3 ΔT(X) as function of speed Ns for Dowson’s bearing.  98  50  55  Model, FE, Khonsari, Song/Hashimoto Temp Dif f Comparisons: Ferron Brg 350  300  ΔTmeasure, ns=67rps =18C  250  ΔT, C  200  150 DT_model DT_FE DT_Khonsari DT_Song/Hashimoto  100  50  0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 Speed, rev/s  (a) Model, FE, Khonsari and Song/Hashimoto ΔT Comparisons. Ferron Brg 60 55 50  Exp. Point DT_model DT_FE DT_Khonsari DT_Song/Hashimoto  45 40  ΔT, C  35  ΔTmeasured, ns=67rps =18  30 25 20 15 10 5 0 10  20  30  40  50  60  Speed, rev/s  (b) Figure 5.4 ΔT(X) as a function of speed Ns for Ferron’s bearing.  99  70  Delta T Comparisons, Costa 600 550 500 450 400  ΔT, C  350 DTModel  300  DT_FE DT_Khonsari  250  DT_Song/Hashimoto  200 150 100 50 0 0  10  20  30  40  50  60  70  80  90 100 110 120 130 140 150 160 170 180 190 200 210  Speed, rev/s  (a) Leakage flow comparisons, Costa 170E-6 160E-6 Measured Leakage Data QL ns=33.3 = 45.8x10-6 m3/s QL ns=50 = 54.2x10-6 m3/s QL ns=66.7 = 62x10-6 m3/s  150E-6 140E-6 130E-6 120E-6 110E-6  QL, m3/s  100E-6 90E-6 80E-6 Exp. Point  70E-6  QL_Model  60E-6  QL_Khonsari QL_Song/Hashimoto  50E-6  QL_FE  40E-6 30E-6 20E-6 10E-6 000E+0 0  10  20  30  40  50  60  70  80  90  100 110 120 130 140 150 160 170 180 190 200 210 Speed, rev/s  (b) Figure 5.5 ΔT(X), QL as function of speed Ns, Costa.  100  Delta T Comparisons, N s =25 rev/s, Dowson's data 30 Measured Data Clearance C = 63.8  25  Exp. Point FE Model Khonsari  ΔT, C  20  15  10  5  0 40E-6  60E-6  80E-6  100E-6 Clearance, C (m)  120E-6  140E-6  (a) Delta T vs Diameter for Dowson's Brg, C=60μm, L=0.0765m 160 150 140 130 120 110  ΔT, C  100 90  ΔT Data (Celsius) Measured = 13.7 Model = 13 Khonsari = 15.8 Song = 39.2 Exp. Point FE Model Khonsari Song/Hashimoto  80 70 60 50 40 30 20 10 0  0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19  Diameter D, (m)  (b) Figure 5.6 (a) ΔT(X) as function of clearance ratio C, Dowson data, and (b) as function of diameter D  101  5.3.3  Leakage Flow QL Predictions  Leakage results were also analyzed for the examples above and presented in Figures 5.7, 5.8 and 5.9. The finite element predictions are closest to the measured values at 25 rev/sec speed. Figure 5.8 presents leakage flow data from Ferron’s bearing (measured at 66.7 rev/s); similar results are evident. Standard error bars are included for the model prediction as well as a select few (to reduce clutter). Model predictions are good in all cases. Hashimoto’s model is significantly off for the range of bearings tested since it is derived for high-speed, short bearings. Another shortcoming with this model ( Q L =  π 4  N s CD 2 ε ) is that it does not specify supply  pressure ps explicitly. Leakage Flow Rate Predictions: Dowson Bearing 150E-6 140E-6 130E-6 120E-6 Ps = 280 kPa QL measured ns=25 = 26.2x10-6 m3/s QL measured ns=33.3 = 30.1x10-6 m3/s  110E-6 100E-6  QL, m 3/s  90E-6 80E-6 70E-6 60E-6  Exp. Points Q_FE Q_Khonsari Q_Song/Hashimoto Q_model  50E-6 40E-6 30E-6 Song/Hashimoto  20E-6 10E-6 000E+0 10  20  30  40  50  60  70  80  90 100 110 120 130 140 150 160 170 180 190 200 210  Speed, rev/s  Figure 5.7 Side leakage predictions for Dowson’s bearing.  102  Leakage Comparisons: Ferron Brg 320E-6 300E-6 280E-6 260E-6  Exp. Point QL_FE QL_Khonsari QL_Song/Hashimoto QL_Model  240E-6 220E-6 200E-6  QL, m 3/s  180E-6 160E-6 140E-6 120E-6 100E-6  QL exp. = 130.6x10 -6 m 3/s  80E-6 60E-6 40E-6 20E-6 000E+0 10  20  30  40  50  60  70  80  90  100  110  120  130  140  150  160  170  180  190  Speed, rev/s  Figure 5.8 Side leakage predictions – Ferron’s bearing. Leakage flow comparisons, Costa 170E-6 160E-6 Measured Leakage Data QL ns=33.3 = 45.8x10-6 m3/s QL ns=50 = 54.2x10-6 m3/s QL ns=66.7 = 62x10-6 m3/s  150E-6 140E-6 130E-6 120E-6 110E-6  QL, m3 /s  100E-6 90E-6 80E-6 Exp. Point  70E-6  QL_Model  60E-6  QL_Khonsari QL_Song/Hashimoto  50E-6  QL_FE  40E-6 30E-6 20E-6 10E-6 000E+0 0  10  20  30  40  50  60  70  80  90  100 110 120 130 140 150 160 170 180 190 200 210  Speed, rev/s  Figure 5.9 Side leakage predictions – Costa’s bearing.  103  200  210  A summary of the data is presented in Table 5.5. The figures in brackets show percentage difference to measured values. The proposed model gives consistently better results than Khonsari and Song models. Table 5.5 Summary of ΔT(X) and QL data comparisons. Note that all leakage flow data are × 10 −6 m3/s, ΔT values are in °C and based on (Tmax-Ti). ΔT(X)  Finite Element  Dowson  13.3 (2.3)  14.2 (-3.4) 16.9 (6) 18.1 23.1 (9.5) (0.4)  43.4 (5.2)  26.2 130 (0.46) 58.9 71.7 (8.7) (15.6)  Ferron Costa QL Dowson Ferron Costa  5.3.4  Model  Khonsari  Song  16 (8.8)  19.1 (29.9)  70.7 (380)  16.0 (23)  19.6 (8.8) 22.2 (11)  21.8 (21) 28.0 (40)  42.2 (133) 122 (500)  41.4 (9.6)  22.5 (14) 134.9 (3.4) 56.9 69.7 (5) (12.4)  Measured  14.7  27.6 (20)  20.0 (54)  35.0 (52)  17.5 (-33) 106.1 (-18.8) 32.6 44.6 54.7  66 (400)  1.05 (-95.9) 5.47 (-96) 1.55 1.61  190 (700)  13  18 20  23  26.2  1.67  45.8  130.6 54.2  62  Variation of C on Performance Parameters ΔT and QL  The model is run to obtain temperature increase ΔT(X) as a function of clearance ratio C in the range 40≤C≤150 μm for Dowson’s bearing. Figure 5.10 shows values at Ns=201 rev/s. The model, FE, and Khonsari predictions follow a similar pattern while Song’s model predictions are well off. Results for side leakage as a function of C for all four models at a speed of 25, 101 and 201 rev/s are shown in figure 5.11. It appears that side leakage is linear at lower speeds and changes to a power function at higher speeds as shown in figure 5.12. A least squares curve fitting was performed in MS Excel with the results QL ∝ C  Ns < 50 rev/s  104  QL ∝ C1.5  50 ≤ Ns ≤ 120 rev/s  QL ∝ C1.3  Ns > 120 rev/s  (5.15)  Included in figure 5.12 are r2 values for each trendline. It has to be mentioned that the bearings analyzed have a diameter close to 100 millimeters. A wider range of sizes needs to be tested to give the results wider applicability.  Leakage Flow at two speeds, Dowson 400E-6 350E-6 300E-6  Leakage, m3/s  250E-6  FE:ns=25 Model:ns=25 Khonsri:ns=25 Song:ns=25 FE:ns=201 Model:ns=201 Khonsari:ns=201 Song:ns=201  200E-6 150E-6 100E-6 50E-6 000E+0 40E-6  60E-6  80E-6 100E-6 Clearance ratio C, (m)  120E-6  Figure 5.10 Side leakage as a function of clearance C, speed 25 to 201 rev/s.  105  140E-6  Delta T Comparisons for Dowson's bearing, N s =201 1200  1000  ΔT, C  800 DTFE:ns=201 DTModel:ns=201 DTKhonsari DTSong/Hashimoto  600  400  200  0 40E-6  60E-6  80E-6 100E-6 Clearance, C (m)  120E-6  140E-6  160E-6  Figure 5.11 ΔT(X) versus clearance ratio C, Dowson’s bearing, Ns=201 rev/s. QL v C, λ=0.75, μi =0.03 Pa.s, Dowson 400.0E-6 350.0E-6  Side Leakage QL, m3/s  300.0E-6  "n_s=25" n_s=101 n_s=200 Linear ("n_s=25") Power (n_s=101) Power (n_s=200)  R² = 0.998  QL α C1.3  250.0E-6 R² = 0.999 200.0E-6 QL α C1.5 150.0E-6 R² = 0.993  100.0E-6  QL α C  50.0E-6 000.0E+0 40.0E-6  60.0E-6  80.0E-6 100.0E-6 120.0E-6 Radial Clearance C, (m)  140.0E-6  Figure 5.12 Side leakage as a function of clearance at three different speeds.  106  160.0E-6  5.3.5  Temperature Increase ΔT, Side Leakage QL as Function of Diameter  The effect of varying diameter D while holding the length L constant for Dowson’s bearing is analyzed. Figure 5.13 shows the temperature increase ΔT(X) while figure 5.14 shows leakage variation. Figure 5.13 is plotted by holding the length of Dowson’s bearing constant and varying the diameter. Song’s model is included for completeness of the discussion although it is clearly unsuitable for the range of bearings tested. Their model is based on the friction force which is calculated by considering the local Reynolds number (Re=ρCU/μ). There appears to be no method of updating the viscosity for use in computation of the Reynolds number in their model. The other models show steady rise with increasing diameter although the gradient ∂ΔT/∂D increases in the region 0.1<D<1.25 after which it decreases.  Delta T vs D, N s=101rps, Dowson 450 400 350  FE@101 Model@101  300  Khonsari@101 Song/Hashi@101  ΔT, C  250 200 150 100 50 0 0.04  0.06  0.08  0.10  0.12  0.14  0.16  Diameter D, (m)  Figure 5.13 ΔT(X) as a function of diameter D.  107  0.18  0.20  Leakage Flow Vs Diameter, Dowson's Data, C=80 mm 7.00E-05  6.00E-05  5.00E-05  Leakage QL, m 3/s  4.00E-05  3.00E-05  Model Khonsari FE Song/Hashimoto  2.00E-05  1.00E-05  0.00E+00 0.04  0.06  0.08  0.10  0.12  0.14  0.16  0.18  0.20  Diameter D, (m)  Figure 5.14 QL as a function of diameter D, Dowson’s bearing, Ns=25 rev/s. The relation is linear as expected at this speed.  5.3.6  Simulation Results with ΔT (f1) and QL (f2) as Objective Functions  The set of equations developed are used to minimize a multi-objective function of temperature increase ΔT(X) (F1) and side leakage QL (F2) using the design variable set  X T = (C, λ, μ ) , (same set used by Hashimoto [34] and Song et al. [36] ). Hashimoto’s model is for high speed (up to 240 rev/s) and short bearings (0.2<λ<0.6) whereas the model developed here is for medium speed (up to 200 rev/s) and medium range (0.4<λ<1.0). The two objective functions are combined using scaling and weighting factors to form one multi-objective function, F(X)=α1β1ΔT+α2β2QL which is minimized subject to linear and non-linear constraints given in equation 5.14(b) and (c). Additional  108  input parameters for Dowson’s bearing to run the optimization model are given in table D.5. Results will be compared to those obtained by Song et al. for high-speed, short bearings where a similar objective function was used (included for comparison only as input parameters are different). Results from both the classical (SQP) and a hybrid method are compared and presented in figures 5.15, 5.16, 5.17 and 5.18. At least three runs with many sub-iterations were performed for each method. Figure 5.15 compares results of the objective function from the two methods at load 10 kN. The hybrid algorithm needed at most two runs to produce a consistent result. SPQ however required at least four runs while adjusting the initial guess to get reasonable predictions. Comparison with Song’s results show a similar pattern to the objective function although Song’s values are lower due to different input parameters and range (see table D.5). Clearance ratio predictions from the model (figure 5.16(b)) show a rapid initial increase and leveling off at speeds above 80 rev/s. This is different from Song’s model prediction (figure 5.16(a)), which shows an initial decline, leveling off, and gradual increase at speeds greater than 120 rev/s. Actual values of C in Song’s model are much lower than model predictions. Length to diameter λ predictions are shown in figure 5.17(b) for the model and 5.17(a) for Song’s model. Both show a decline to the lower bound of λ, the model showing a faster descend than Song’s model. Thus a similar pattern is exhibited by both models. Finally viscosity variation μ for both models show predictions in the lower bound values as shown in figure 5.18 (a) and (b). In summary, the following points are evident: •  Objective function predictions F(X) show a smooth increase for both models  109  •  Clearance ratio requirements for wider bearings appear to be higher than shorter bearings  •  Length to diameter ratio λ decreases to the lower bound in both cases, declining faster for wider bearings.  •  Viscosity remains at the lower bound in both cases. Obviously viscosity values should be as low as possible to minimize viscous dissipation.  The next section extends the design variable set to four ( X T = (C , λ , μ , p s ) ) by including supply pressure variation; the optimization work primarily focuses on pressure fed bearings. The multi-objective function is also reformulated by considering minimization of the power loss PL and leakage QL ; these two are axiomatically independent [39].  110  Objective Function: SQP vs GA Objective Funcition f(X)  18 16 14 12 10  SQP  8  GeneticAlg  6 4 2 50  100  150  200  250  Speed, rev/s  (a) Objective Function f(X): SQP vs Hybrid 60  50  W = 10 kN  f(X)  40  30  f(X): Hybrid SQP: 7 Attempts  20  10  0 0  10  20  30  40  50  60  70  80  90  100 110 120 130 140 150 160 170 180 190 200  Speed, rps  Figure 5.15 (a) Objective function from Song et al. [36], and (b) F(X) predictions from SQP and hybrid algorithm, Wd=10kN for both cases.  111  Radial Clearance, μm  Radial Clearance C, μm  100 90 80 70  SQP GA  60 50 40 50  100  150  200  250  Speed, rev/s  (a) Radial Clearance C : SQP vs Hybrid 250E-6 W = 10 kN  Radial Clearance C, m  200E-6 C: SQP C: Hybrid 150E-6  100E-6  50E-6  000E+0 0  10  20  30  40  50  60  70  80  90  100 110 120 130 140 150 160 170 180 190 200 210  Speed, rps  (b) Figure 5.16 (a) Clearance ratio C as a function of speed for Song’s model (b) Model predictions, Wd=20 kN  112  Length to Diameter Variation λ  Length to Diameter λ  0.30  0.25  0.20  SQP  0.15  GA 0.10 50  100  150 Speed, rev/s  200  250  (a) Length to Diameter Ratio λ Comparisons 1.2  W = 10 kN 1  λ  0.8  0.6  0.4  1:SQP: 7 attemps 0.2  2: Hybrid  0 0  10  20  30  40  50  60  70  80  90  100  110  120  130  140  150  160  170  180  190  200  Speed, rps  Figure 5.17 (a) Length to diameter ratio λ as function of speed, Song’s model, and (b) Length to diameter ratio λ as function of speed for SQP and hybrid algorithm for developed model.  113  Viscosity μ, Pa.s  Viscosity Variation μ 1.10 1.08 1.06 1.04 1.02 1.00 0.98 0.96 0.94 0.92 0.90  SQP GA  50  100  150  200  250  Speed, rev/s  (a) Viscosity μ Comparisons 0.033  0.03 W = 10 kN 0.027 2: Hybrid SQP: 7 Attempts  Visocosity μ, (Pa.s)  0.024  0.021  0.018  0.015  0.012  0.009  0.006 0  10  20  30  40  50  60  70  80  90  100 110 120 130 140 150 160 170 180 190 200  Speed, rps  (b) Figure 5.18(a) Viscosity variation, Song/Hashimoto’s model [36] (b) Viscosity variation versus speed, hybrid and SQP.  114  5.4  Power Loss PL (F1) and Leakage QL (F2) as Objective Functions  5.4.1  Power Loss Predictions from Model  Power loss comparisons for the proposed model, Jang and Khonsari’s [55], finite element, and Hashimoto/Song’s models [34,36] is shown in figure 5.19 for Dowson’s bearing [59]. Five measured data points and standard error bars are included; these show the consistently better prediction from the model. Figure 5.20 presents simulation results for Costa’s data [69]. Both figures show the relative accuracy of the proposed model for the range of bearings tested (0.4<λ<1.5 and 0<Ns<200 rev/s). Model predictions are well within the 95% confidence level. Another advantage is the ease with which the model can be implemented. Equation set 5.9 computes the coefficient of friction f, friction force F, and power loss PL. Power Loss Comparisons - Dowson's Data 1800 1700 1600 Exp. Data Speed, rps Loss, W 11.7 13.3 16.7 20.8 25  1500 1400 1300 1200  Power Loss, W  1100 1000  Power 153.1 189.2 280.8 406.4 532  900 800 700 600 500 400  Exp. Points FE Model Khonsari Song/Hashimoto  300 200 100 0 5  7  9  11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51  Speed, rev/s  Figure 5.19 Power loss prediction – Dowson’s data.  115  Power Loss Comparisons, Costa Bearing, 70 kPa 9000 8500 8000  Measured Data PowerLoss ns=33rps =612 W PowerLoss ns=50rps =1125 W PowerLoss ns=67rps =1663 W  7500 7000 6500  Power Loss, W  6000 5500 5000 4500 4000 3500 3000 2500 Exp. Point FE Model Khonsari Song/Hashimoto  2000 1500 1000 500 0 0  10  20  30  40  50  60  70  80  90  100 110 120 130 140 150 160 170 180 190 200 210  Speed, rev/s  Figure 5.20 Power loss predictions – Costa’s data.  5.4.2  Determination of weighting parameters and scaling factors  The method of weighting/scaling factors is used to determine optimum values for the design variables. Input data shown in Tables D.1 and D.5 is used. Values of leakage QL were weighted appropriately by using a suitable value for β2 while the scaling factors were selected after comparing the absolute values of the objective function F(X) as a function of the ratio α1/α2. Three values of the ratio were tried, namely 5, 1, and 1/5. It is assumed the objective functions f1 and f2 are of equal importance and this means the magnitude of the two should be equal or nearly so. The effect of weighting factor β2 of 105 is to transform leakage values, which are in the order 10-5 m3/s to the same order of magnitude as power loss. Values of f1 and f2 are adjusted to be of equal magnitude by  116  selecting α1/α2=1/5, β1=1/400, β2=105 as shown in figure 5.22a. It is important that values of the objective functions be similar if they are of equal importance (somewhat like ‘comparing oranges with oranges’). It must be emphasized that the ultimate choice of the weighting parameters depends on the relative importance of each objective function f1 and f2 [34]. These graphs show that power loss f1 increases more rapidly than leakage f2 thus having a greater influence on the objective function value F(X). Depending on the relative importance of minimizing the power loss or leakage, figure 5.21 gives an indication on what route to take. An assumption is made here that the goal is to obtain the lowest possible value of the objective function and use α1/α2=1/5. The effect of scaling factors β1 and β2 with α1/α2=1/5 is shown in figures 5.22, 5.23 and 5.24. Also included are the values of F1, F2 and the design variable sets (C,λ,μ,ps). One observation is increasing clearance ratio reduces power loss but increases side leakage significantly. Figure 5.21 shows results with β1=1/400 and β2=105 which has the effect reducing power loss contribution to F(X). The net result is minimizing side leakage has more weight and clearance ratio is at the minimum value in 5.21c while length to diameter ratio λ is at the upper limit of 1. Supply pressure predictions show rapid increase with increasing speed to the upper limit at a speed of 40 rev/s. Viscosity is virtually at the minimum allowable in all cases.  The trends change significantly  however, when β1 is increased to 1/35 in figure 5.22, making minimization of power loss (F1) more important (figure 5.22a). Supply pressure increases rapidly to the upper limit as before, however higher clearance ratio values are predicted. Length to diameter values shift lower while viscosity remains at the lower values. Significant changes occur with  β1=1 (figure 5.23) as supply pressure increases less rapidly, viscosity remains at the  117  lower limit, clearance ratio is initially high and reduces, while length to diameter ratio λ reduces further to values in the region of 0.65 at higher speeds.  Objective Function Values for Different Scaling Factors  3000 α1/α2=5 α1/α2=1 α1/α2=1/5  Obj. Function f(x)  2500  2000  1500  1000  500  0 0  20  40  60  80  100  120  140  160  180  Speed, rps  Figure 5.21 f(X) as function of α1/α2, scaling parameters: β1=1, β2=105  118  200  Objective function F(X), F 1 , F 2 for β 1 =1/400 Objective function  30 f(X) f_1 f_2  25 20 15 10 5 0 0  20  40  60  80 100 120 140 160 180 200  Speed, rps  (a) Supply pressure ps variation 400E+3 350E+3 300E+3 250E+3 200E+3 150E+3 100E+3 50E+3 000E+0  Variation of clearance ratio C Radial Clearance C, m  42E-6  ps, kPa  p_s  0  40E-6 38E-6  34E-6 32E-6 30E-6  20 40 60 80 100 120 140 160 180 200  0  Speed, rps  20  40  60  80 100 120 140 160 180 200  Speed, rps  (b)  (c)  Length to diameter ratio λ  Viscosity variation μ  1.1  0.035 0.030 0.025  μ  0.020  λ  Viscosity μ, (Pa.s)  C  36E-6  0.015  1  0.010  λ  0.005 0.000 0  20 40 60 80 100 120 140 160 180 200  0.9 0  Speed, rps  20 40 60 80 100 120 140 160 180 200  Speed, rps  (d) (e) Figure 5.22 Model predictions with α1/α2=1/5, β1=1/400, β2=105  119  Objective function F(X) , F 1 and F 2 140  Objective function f  120 100 80  f(X)  60  f_1  40  f_2  20 0 0  20  40  60  80  100 120 140 160  180 200  Speed, rps  (a) Variation of clearance ratio C  400  48E-6  350  47E-6  Radial Clearance C, m  Supply press. ps , kPa  Variation of supply pressure p s  300 250 200 p_s  150 100 50  46E-6 45E-6 44E-6 43E-6  C  42E-6 41E-6 40E-6 39E-6  0 0  20 40  0  60 80 100 120 140 160 180 200 Speed, rps  20 40 60 80 100 120 140 160 180 200 Speed, rps  (c)  (b)  Variation of λ: Hybrid  Variation of inlet viscosity μ  1.05 1  0.01 0.008  μ  λ  0.95  λ  Viscosity μ, (Pa.s)  0.012  0.006  0.9  0.004  0.85  0.002  0.8 0  0 0  20 40 60 80 100 120 140 160 180 200  20 40 60 80 100 120 140 160 180 200  Speed, rps  (e)  Speed, rps  (d) Figure 5.23 Model predictions with α1/α2=1/5, β1=1/35, β2=105  120  Obj. function  Objective functions F(X), F 1 , F 2  4500 4000 3500 3000 2500 2000 1500 1000 500 0  f(X) f_1 f_2  0  20  40  60  80  100 120 140 160 180 200  Speed, rps  (a)  Variation of clearance ratio C 8.5E-05 8.0E-05 7.5E-05 7.0E-05 6.5E-05  C, m  ps, kPa  Variation of supply pressure ps 400 350 300 250 200 150 100 50 0  p_s  6.0E-05 5.5E-05  C  5.0E-05 4.5E-05 4.0E-05  0  20 40 60 80 100 120 140 160 180 200  0  Speed, rps  Speed, rps  (b)  (c) Length to diameter ratio λ  Variation of inlet viscosity μ  1  0.012  0.9  0.01 0.008  0.8  μ  λ  Viscosity μ, (Pa.s)  20 40 60 80 100 120 140 160 180 200  λ  0.7  0.006 0.6  0.004  0.5  0.002  0.4 0  0 0  20 40 60 80 100 120 140 160 180 200  20  40 60 80 100 120 140 160 180 200  Speed, rps  Speed, rps  (e)  (d)  Figure 5.24 Model predictions with α1/α2=1/5, β1=1, β2=105  121  5.4.3  Pareto Optimal Front Method  The method of Pareto diagrams eliminates the need for scaling and weighting factors [40,42]. Pareto fronts for a load of W=10 kN are developed. The two objective functions are plotted on the axes and a series of non-inferior points plotted. These charts can be very effective in helping designers choose design variable sets depending on the relative importance they place on each objective function. Pareto optimal fronts at common speeds for these bearings (1500 rev/min (25 rev/s), 2000 (33.3), 3000 (50) and 4000 (66.7)) are shown in figures 5.25 and 5.26. Included in each chart is a boxed value for the recommended operating point. The design variable set for this point is also shown. The power loss range is small in these plots because the Pareto fronts are developed by changing the supply pressure and determine its effect on both objective functions. The Pareto optimal fronts all satisfy the condition in equation 5.12. The Pareto table at 4000 rpm is shown in Table 5.6.  122  26.70E-6  Pareto Optimal Front, N s =25 rps  26.60E-6  Leakage m3/s  26.50E-6 26.40E-6  {6.94E-5,0.734,0.01,75E3}  26.30E-6 26.20E-6 26.10E-6 26.00E-6 25.90E-6 25.80E-6 25.70E-6 418.5  419  419.5  420  420.5  Power Loss, W  (a) Pareto Optimal Front, Ns=33.3 rps 31.80E-6  Leakage, m3/s  31.60E-6  31.40E-6  31.20E-6 {65.2E-6,0.697,0.01,107E3} 31.00E-6  30.80E-6  30.60E-6 583.6  584.1  584.6  Power Loss, W  (b) Figure 5.25 Pareto optimal front for speeds of (a)1500 rpm and (b) 2000 rpm.  123  Pareto Optimal Front, N=50 rps 42.40E-6  Leakage, m 3/s  42.30E-6 42.20E-6 42.10E-6  {6.11E-5,0.664,0.01,184E3}  42.00E-6 41.90E-6 41.80E-6 931.7  931.85  932  932.15  932.3  932.45  Power Loss, W  (a) Pareto Optimal Front, Ns=66.7 rps 54.20E-6 54.10E-6 54.00E-6  Leakage, m 3/s  53.90E-6 53.80E-6 53.70E-6 53.60E-6 {59.9E-6,0.652,0.01,250E3}  53.50E-6 53.40E-6 53.30E-6 53.20E-6 1257.3  1257.45  1257.6  1257.75  1257.9  1258.05  Power Loss, W  (b) Figure 5.26 Pareto optimal front for (a) 3000 rpm and (b) 4000 rpm. Power loss in Watts.  124  Table 5.6 Typical Pareto Table for Bearing at 4000 rpm, W=10kN Length to Power Loss  Leakage  Clearance  Diameter  Supply  PL, Watts  QL (m3/s)  C (μm)  ratio λ  pressure ps  1257.354 1257.354 1257.381 1257.419 1257.420 1257.470 1257.540 1257.560 1257.607 1257.652 1257.687 1257.721 1257.776 1257.815 1257.877 1257.916  54.138E-6 54.051E-6 53.950E-6 53.801E-6 53.796E-6 53.670E-6 53.546E-6 53.516E-6 53.464E-6 53.421E-6 53.394E-6 53.371E-6 53.339E-6 53.323E-6 53.298E-6 53.288E-6  60.050 60.051 60.031 60.002 60.002 59.965 59.913 59.898 59.864 59.830 59.805 59.780 59.740 59.711 59.666 59.638  0.65197 0.65205 0.65201 0.65203 0.65203 0.65200 0.65195 0.65192 0.65188 0.65182 0.65178 0.65173 0.65166 0.65162 0.65151 0.65146  140.0E+3 150.0E+3 160.0E+3 179.1E+3 180.0E+3 200.0E+3 225.0E+3 235.0E+3 250.0E+3 260.0E+3 270.0E+3 279.5E+3 293.9E+3 305.0E+3 320.0E+3 330.0E+3  5.5  Extended Optimization Model  5.5.1  Leakage Flow Equation and Validation  The flow equation presented by Zengeya and Gadala [76] and detailed in section 5.3 assumes flow is a function of aspect ratio L/D and eccentricity ratio ε, (QL=f(λ,ε)) only; it does not explicitly consider the effect of supply pressure and oil groove geometry on leakage flow rate. The method proposed by Martin [70] considers leakage flow as the sum of a hydrodynamic pressure-induced and velocity-induced components. The drawback of Martin’s work is in the many considerations of bearing inlet conditions (oil hole/oil groove, film thickness at inlet port hg, and parameter fg which is a function of the oil groove or hole dimensions). This justifies the need for formulating a simplified leakage flow model with less complex considerations while maintaining accurate predictions. The  125  work presented by Zengeya and Gadala [76] is used as a base and expanded here to take the effect of groove geometry Lg/L, groove location θ, groove length to diameter ratio Lg Lg ⎞ ⎛ Lg/D into account, i.e., Q L = F ⎜⎜ ε , λ , p s , θ , , ⎟⎟ . Based on extensive numerical D L ⎠ ⎝  simulations using the streamline upwind Petrov-Galerkin finite element method, the following equation for leakage flow rate is proposed: Lg ⎛ Q L = (1 + 1.6 p s )(1 + 0.17θ )⎜⎜1 + 0.3 D ⎝  L ⎞⎛ ⎟⎟⎜⎜1.15 g L ⎠⎝  ⎞ ⎟⎟ exp 0.2 − 0.5λ1.5 + 1.05 ln ε ⎠  [ (  ps where p s , the non-dimensional supply pressure p s = 2πμN s  )]  (5.16)  2  ⎛C ⎞ ⎜ ⎟ is determined using ⎝R⎠  the fluid pressure ps, the fluid viscosity μ, the journal speed Ns (rev/s), and the radial clearance C. Other parameters in equation 5.16 include the length to diameter ratio λ = L/D, the groove length Lg, and the groove location θ. It should be noted that θ can be  positive or negative as shown in figure 4.1a; it varies between ±π/2. The exponential part of the equation is adopted from the work of reference 76. Depending on the nature of the bearing set-up, appropriate terms in equation 5.16 may be ignored. Justification for the choice/form of p s and θ (linear) in equations 5.16 follows. Figures 5.27 to 5.31 show the linear relationship between leakage flow and supply pressure, groove location and eccentricity ratio for some bearing data. Figure 5.27a illustrates the effect of supply pressure on flow rate for Claro’s data [79] at eccentricity ε=0.7, speed Ns=8.6 rev/s and supply pressure variation from 8 to 60 kPa. Maximum error of 17 per  cent from experimental data is better than 25 per cent reported in Claro’s work (error=(predicted-experimental)/experimental value)x100). Standard error bars also  126  reveal two data points outside the 95% confidence level. Figure 5.27b shows model predictions on Costa’s data [69] at speed Ns=50 rev/s, eccentricity ε=0.5 and supply pressure values ranging from 70 to 300 kPa. Maximum error is 19 per cent with two points outside the 95% confidence level. Figure 5.28 presents flow as a function of speed at two supply pressures for Dowson’s data [80,81]. Predictions from the model show a maximum error of 14 per cent while Dowson’s is 21. Both models fall outside the 95% confidence level for the experimental data at ps=170kPa. Figure 5.29 shows the effect of both supply pressure and eccentricity ratio on leakage flow for Dowson [80] whose model predictions are included for comparison. Both models are accurate at zero pressure; however the model developed here shows better prediction of the data as indicated by determining the variance of each in table 5.7. The goodness of fit, however, is poor at  1.0. The accuracy of the model at two supply pressures is tested in figure  5.30 while figure 5.31 presents non-dimensional leakage flow as a function of eccentricity ratio from Boncompain [15] and Ferron [75]. Two other predictions are included for comparison. Error bars in figure 31a indicate all models fail to accurately predict higher eccentricity values.  127  Measured vs Predicted Leakage, Claro et al., θ=+90o 0.6  0.5  Non-dim. leakage  QL 0.4 Measured  0.3  Model  0.2  0.1  0 0  10  20  30  40  50  60  70  Supply pressure, kPa  (a) Predicted vs Measure values, Costa (p s=0.03,0.04,0.09,0.13,0.19 kPa, θ=+30)  0.7 0.6  Non-dim. leakage  QL  0.5 0.4  Model Measured  0.3 0.2 0.1 0 0  0.05  0.1  0.15  0.2  Non-dim supply presure p s  (b) Figure 5.27 (a) Leakage flow prediction and data as a function of supply pressure for Claro’s data [79]. Groove at +90° to the load line. Other bearing parameters as in Table D2. (b) Data from Costa et al. [69]. Groove location: +30°  128  Predicted, Measured, and Dowson Leakage Data 2.50E-05  Leakage flow QL, m 3/s  2.00E-05  1.50E-05  Measured:ps=170kPa Q_LDowson:ps=170kPa QMeasured:ps=0 Q_LModel:ps=0 Q_LDowson:ps=0 Q_LModel:ps=170kPa Q_LKhonsari:ps=0 Q_LKhonsari:ps=170kPa  1.00E-05  5.00E-06  0.00E+00  0  5  10 Speed, rev/s  15  20  Figure 5.28 Leakage flow as function of supply pressure, Costa [69]. Ns=50 rpm, ε=0.5. Non-dimensional Flow: Predicted vs Measured Values 1.40 Measured:ps=0 Qmodel:ps=0  1.20  Measured@ps=0.5 Qmodel:ps=0.5 Mesured@ps=1  Non-dimensional flow  1.00  Qmodel:ps=1 Dowson:ps=0 Dowson:ps=0.5  0.80  Dowson:ps=1  0.60  0.40  0.20  0.00 0.15  0.25  0.35  0.45 0.55 Eccentricity ratio ε  0.65  0.75  0.85  Figure 5.29 Non-dimensional leakage flow as a function of eccentricity ratio for Boncompain [15] and Ferron [75]. 129  Table 5.7 Comparison of variance between model and Dowson’s predictions.  Variance  ps  Model Dowson 0 0.0022 0.0015 0.5 0.0065 0.0157 1.0 0.0121 0.0176 Leakage Flow Vs Speed, Dowson/Taylor 6.0E-05 Model predictions (measured data in brackets) -6 -6 3 p s =0: Q L,ns=13.3=5.4x10 (5.2x10 ) m /s  5.0E-05  -6  -6  3  3  Leakage flow, m /s  p s =170kPa: Q L,ns=13.3=22.1x10 (21x10 ) m /s  4.0E-05 3.0E-05 ps170  2.0E-05  ps=0  1.0E-05 0.0E+00 0  20  40  60  80  100  120  140  160  180  200  Speed, rps  Figure 5.30 Leakage flow as a function of speed ps=0 and 170 kPa, Dowson’s [80].  130  Leakage Flow vs Eccentricty, 1500 rpm  Leakage Flow vs Eccentricity, 2000 rpm  1.2 Measured  1  Khonsari  0.8  model khonsari Martin Measured  QL 0.8  Martin  Non-dimensional f low  Non-dimensional flow  QL  1.0  Model  0.6 0.4 0.2  0.6 0.4 0.2 0.0  0 0.30  0.40  0.50  0.60 0.70 Eccentricity ε  0.80  0.15  0.90  0.25  0.35  0.45  0.55  0.65  0.75  0.85  Eccentricity ε  (b)  (a)  Leakage Flow vs Eccentricity, 3000 rpm  Leakage Flow vs Eccentricity, 4000 rpm 0.8  0.9 model  QL0.8  QL  khonsari  0.6  martin  Non-dimensional f low  0.7  Non-dimensional flow  Exp. Data  0.5 0.4 0.3 0.2  0.7 0.6 0.5  model khonsari martin Exp. data  0.4 0.3 0.2 0.1  0.1 0.0  0 0.30  0.40  0.50  0.60  0.70  0.80  Eccentricity ε  0.20  0.30  0.40  0.50  0.60  0.70  Eccentricity ε  (d)  (c)  Figure 5.31 Non-dimensional leakage flow QL as a function of eccentricity ratio for Boncompain [15] and 75 [75].  5.5.2  Power Loss Equation and Validation  Power loss determination follows a similar procedure presented in section 5.4. The coefficient of friction parameter (R/C)f in reference 76 is expanded to include effects of  131  supply pressure and groove position θ  and with the aid of extensive numerical  simulations and curve-fitting techniques: ⎡ ⎛ ⎞⎤ R 2.73622 − 4.97838 ε ⎟⎟⎥ f = 1 − 0.25θ 2 + 0.05θ (1 + 0.6 p s )⎢exp⎜⎜1.9994 + C λ ⎠⎦ ⎣ ⎝  (  )  (5.17)  where f is the coefficient of friction and ε is the eccentricity ratio, ε = e/C. The above equation introduces the effects of supply pressure p s and groove position θ in an explicit way. It will become evident that the groove position plays a significant role in calculating the friction and the subsequent torque and power values from validation simulations. The effect of supply pressure, although more significant at low eccentricity ratios [80], is included in both equations 5.16 and 5.17. The quadratic nature of coefficient of friction as a function of groove location in equation 5.17 is evident in figure 5.32 for Costa’s data speed Ns=50 rev/s, eccentricity ε=0.21 and 0.5 (load W is 2 and 8 kN respectively). The magnitude of each coefficient is determined empirically and averaged over the bearing data. Maximum error is 9.4 per cent. Variation of torque as a function the non-dimensional supply pressure is illustrated in figure 5.33, Costa [69]. Supply pressure varies from 50 kPa to 300 kPa (non-dimensional values from 0.04 – 0.23) at two different bearing loads (2 kN, ε=0.2 and 8 kN, ε=0.5) which is in the light to medium bearing load range [80]. Maximum error is 14 per cent. Figure 5.34 confirms the accuracy of the model in predicting Costa’s data at three speed values (Ns=33.3, 50 and 66.7 rev/s, θ=+30°) and eccentricity of 0.5. Values in brackets are experimental.  132  Predicted Vs Measured Torque: ps =70kPa 4.50 Data Measured at 2kN: 2.7 N.m Predicted at 2kN: 2.9 N.m (9.4%) Measured at 8kN: 3.7 N.m Predicted at 8kN: 4.0 N.m (6.8%)  TvTheta:Model_W=2kN Measured@2kN TvTheta:Mode_W=8kN  4.00  Torque, N.m  3.50  Measured@8kN  3.00  2.50  -0.6  -0.5  -0.4  -0.3  -0.2  2.00 -0.1 0 0.1 0.2 Groove position, radians  0.3  0.4  0.5  0.6  Figure 5.32 Torque predictions as a function of groove position, Costa [69]. Supply pressure of 70 kPa. Predicted vs measured torque values, θ=-30o 4.00 3.50  Torque, N.m  3.00 2.50 2.00  Model@2kN Measured@2kN Model@8kN Measured@8kN  1.50 1.00 0.50 0.00 0  0.05  0.1  0.15  0.2  0.25  Non-dimensional supply pressure  Figure 5.33 Torque predictions as a function of non-dimensional supply pressure, Costa et al. [69]. 133  Torque T as function of speed, θ=+30o , ps=70kPa  6.00 Measured data (Predictions in brackets): TNs=33.3=3.1 N.m (3.08) TNs=50= 3.6 N.m (3.66) TNs=66.7=4.05 N.m (4.02)  5.00  Exp. Data T  Torque, N.m  4.00  3.00  2.00  1.00  0.00 0  10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 Speed, rev/s  Figure 5.34 Variation of torque as function of speed, Costa [69]. Three data points are shown.  5.5.3  Optimization Simulations  The program structure follows the steps in Appendix C and is implemented in Matlab [51] with input data in tables D5 (Dowson) and D6. The two objective functions are chosen to be: F1 the non-dimensional friction parameter ((R/C)f) and F2 the leakage flow (QL). A weighting factor β2 of 1x105 is used to make the two objective functions of the  same order of magnitude while α1 and α2 are both set to 1. Figure 5.35 shows the variation of design variables X T = (C , λ , μ , p s , θ ) as a function of speed. Optimum radial clearance increases with speed to a peak at approximately 30 rev/s before decreasing and leveling off close to 60 μm. It is interesting to compare this value 134  Optimum L/D λ ratios  Optimum radial clearance C 1.4E-04  1  Radial clearance C  1.2E-04  0.9  1.0E-04  λ  0.8 L/D ratio λ  C  0.7  8.0E-05  0.6  6.0E-05  0.5  4.0E-05  0.4  2.0E-05  0.3 0.0E+00  0 0  20  40  20 40 60 80 100 120 140 160 180 200 Speed, rps  60  80 100 120 140 160 180 200  Speed, rpm  (b)  (a)  Optimum supply pressure p_s  Optimum viscosity μ  60 Supply pressure, kPa  Viscisity, Pa.s  0.015  0.01 μ  0.005  50 40 p_s  30 20 10 0  0 0  20 40 60 80 100 120 140 160 180 200  0  20 40 60 80 100 120 140 160 180 200  Speed, rps  Speed, rps  (c)  (d) Optimum groove position θ 0.0  Groove position θ  -0.2 0  20  40  60  80 100 120 140 160 180 200  -0.4 -0.6 -0.8  θ  -1.0 -1.2 -1.4 -1.6 -1.8  Speed, rps  (e) Figure 5.35 Variation of XT as a function of speed, W=10kN. Model predictions with α1/α2=1, β1=1, β2=105  135  with 138.8µm used by Dowson [80]. The predicted value is obtained with groove position θ=-π/2. Their data shows bearing speeds up to 16.7 rev/s. Another interesting observation is evident when figure 5.35a is compared to figure 11c in reference 76 where the groove is positioned at the load line. The predicted peak clearance of 120 μm at 15 rev/s in figure 5.35a is significantly higher than a peak of 79 μm at 10 rev/s in figure 11c, reference 76. Corresponding C values at each speed are generally higher in figure 5.35a, suggesting optimum groove location allows for higher optimum clearance ratios. Length to diameter ratio values (figure 5.35b) decrease from a high of 0.94 and levels off at 0.67. This pattern is similar for both bearings simulated and is also comparable to figure 11c, reference 76. Optimum viscosity is at the lower bound as expected. This is because lower viscosity lowers power loss as well as increases leakage for the same set of bearing input variables. Similar findings are reported by other researchers [34,36]. Supply pressure is predicted to be at the lower bound for all speeds. This is significantly different from predictions in figure 11b, reference 76 where it increases from the lower to upper bounds. Groove location appears to have a significant effect on supply pressure. Lower values can be used when optimum location is used in journal bearing design. Groove location prediction is interesting in that the model indicates optimum placing of the groove at the lower bound (-π/2), figure 5.35e. Although this position has the least leakage flow, it has significant power loss minimization. Obviously no consideration of whether this position is in the rupture zone and what effect that it has on mechanical stresses is taken. For convenience most journal bearings are designed with groove  136  location in line with the load. The other potential beneficial position, although not shown in the figures, is at the upper bound (+π/2) which gives slightly higher F(X) values. Figure 5.36 shows variation of the design variables as a function of speed for a load of W=700N. The results are similar to those in figure 5.35 although peak C value of 107  μm occurs at higher speed (25 rev/s). Length to diameter ratio levels off at a lower value of 0.45 in figure 5.36b compared to 0.67 in figure 5.35b. Pareto charts present an alternative method of showing optimum design points. A typical pareto optimal front for load W=10kN and speed 50 rev/s is shown in figure 5.37a. In figure 5.37b the load is 700N and speed 25 rev/s. Two finite element predictions with the boxed XT values are included, showing the accuracy of the optimization model. The graph shows significant power loss variation. Pareto charts present a set of nondominated solutions, thus each point on the front has its own merits. A wide range of these charts can be produced at small speed increments. Designers can then make use of these charts in design decisions. An interesting observation is in the power loss values predicted in figure 5.37a (740W) compared to figure 13a (932W), reference 76 for similar input values. Significant power loss savings can be realized with proper groove location.  137  Optim um radial clearance C as function of s peed  Variation of L/D ratio λ  1.2E-04  1.0 0.9 0.8  8.0E-05  0.7  C  λ  0.6 6.0E-05  λ  Radial clearance C, m  1.0E-04  0.5 0.4  4.0E-05  0.3 2.0E-05  0.2 0.1  0.0E+00 0  20  40  60  0.0  80 100 120 140 160 180 200  0  20  40  60  80 100 120 140 160 180 200  Speed, rps  Speed, rps  (a)  (b) Optim um s upply pres s ure ps  0.01  50  Supply pressure. kPa  60  0.008 0.006 μ  0.004 0.002  40  ps 30 20 10  0  0 0  20  40  60  80  100 120 140 160 180 200  0  20  40  60  80  100 120 140 160 180 200  Speed, rps  Speed, rps  (c)  (d) Optimum groove location θ  0.0 -0.2 0  Groove location, radians  Viscosity, Pa.s  Optimum viscosity μ as function of speed 0.012  -0.4  20  40  60  80  100 120 140 160 180 200  Speed, rps  -0.6 -0.8 -1.0  θ  -1.2 -1.4 -1.6 -1.8  (e) Figure 5.36 Variation of XT as a function of speed, load W=700N.  138  Pareto front, W= 10kN, N s =50 rev/s 3.70E-05  3.50E-05  FE  63.3e-6,0.818,01,50e3,-1.571  3  Leakage flow, m /s  Pfront  68.1e-6,0.699,01,50e3,-1.571  3.60E-05  3.40E-05 3.30E-05 3.20E-05 3.10E-05 3.00E-05 710  720  730  740 750 Power Loss, W  760  770  780  (a) Pareto Chart, Speed=1500 rpm, Load=700N 1.2E-05  1.0E-05  Ns=1500rpm  3  Leakage Flow, m /s  FE  8.0E-06  6.0E-06  4.72e-5,0.47,01,50e3,-1.571  4.0E-06  1.16e-4,0.728,01,50e3,-1.571  2.0E-06  0.0E+00 35  45  55  65  75  85  Power Loss, W  (b) Figure 5.37(a) Pareto optimal front for speeds of 1500 rpm, W=10kN; (b) Pareto optimal front, W=700N  139  5.6  Interpretation of Optimization Results  A brief overview of the optimization results is appropriate here. Results from the four variable model are analyzed in more detail and compared to five variable ones.  5.6.1  Pareto Optimal Fronts with X  , , ,  and X  , , ,  ,  The Pareto charts in figure 5.25 show a variation of power PL and leakage flow QL of less than one and four percent respectively. The power loss variation is small, and in engineering terms might be interpreted to mean other factors have more influence in determining bearing performance. Noting that for this set the load is in line with the oil groove, it is clear this position limits the extent to which power loss and side leakage can be influenced. The situation is however significantly different when groove location θ is included in the design set. As shown figure 5.37a power loss variation is 8% while leakage flow varies 23% which gives a significant influence on the two performance parameters. Table 5.7 summarizes the results from the two sets. Further analysis of the results appear in reference 76. Table 5.8 Variation of Performance Parameters PL and QL X X  , , , , , ,    ,  Characteristic Power Loss Leakage flow Power Loss Leakage flow  140  Variation % 0.2 1 8 23  5.6.2  Uncertainty in the Results  The optimization model is based on fundamental engineering principles of lubrication flow and temperature increase as proposed by the respected research teams Jang and Khonsari [55] and Khonsari et al. [74]. There is no reason to doubt the validity of the temperature rise parameters in reference 76 or the form of the flow equation in reference 76. What might be questioned is the accuracy of the data used. Dowson’s data [59] is regarded as a benchmark in lubrication circles while Costa et al. [69] present a wide variety of performance data which is also well taken in the industry. Unfortunately both sets of researchers do not present any error analysis with their data. The finite element simulations presented in earlier chapters (and Appendix E) were shown to be accurate in predicting experimental data from Boncompain [15], Ferron [75], Costa [69] and others [72,73,79,80,81]. Prediction results were compared to available models and shown to be at least as good or better [73].  5.7  Chapter Summary  An Empirically Derived Optimization Model (EDOM) is proposed for length to diameter ratio 0.4≤λ≤1.0. Both minimization of the temperature increase ΔT and leakage QL as well as power loss PL and leakage QL are analyzed. The objective function minimizing power loss PL and side flow QL is preferable as it minimizes two independent functions. The proposed approach can be used by design engineers who need simple and effective models to help them choose optimum design variables using readily available optimization programs such as Matlab’s Optimization Toolbox. Such an approach presents advantages over more complex methods such as Hirani’s model in that no complex integration techniques are required.  141  Model predictions for temperature increase, leakage flow, and power loss are accurate in comparison to models in the literature. The equations are empirically derived and make the need for complex integrations unnecessary. The method of weighting functions/scaling parameters as well as Pareto optimal fronts are used. The benefits of each method are assessed. Pareto optimal diagrams eliminate the need to determine the weights and scaling required in the first method and should be used more often in journal bearing design. The model is extended to include groove position and supply pressure computations explicitly by introducing terms width to diameter ratio (wg/D) and groove length to bearing length ratio (Lg/L). Maximum error of 15% on leakage is observed from the proposed model. The model is comprehensive yet simple to implement, and depending on the nature of the problem some terms in the equations may be ignored. Power loss and torque predictions are within 10% of experimental data. These results show the effectiveness of the model and the ease of use in bearing optimization problems. Model results show that groove location has a significant effect on the optimization of the bearing. The model predicts groove location should be at the lower bound (-π/2). Higher values of clearance ratio C may be used when optimum groove position is chosen. Significant power loss savings may be realized with appropriate location of oil groove.  142  Chapter 6 Summary, Conclusions and Future Work  6.1  Summary and Conclusions  Journal bearing design and optimization requires accurate and reliable prediction of thermal fields due to viscous dissipation. The effect of temperature increase on load carrying capacity of the bearing can be dramatic due to a decrease in the viscosity of the film. Several thermo-hydrodynamic (THD) models in the literature focus on formulating a robust, rapid converging THD model that can be used in optimization where many iterations at varying speeds are required. This requirement presents one of the main challenges for researchers and is the focus of this work. A robust and effective FEA procedure is proposed and a novel and simple model for journal bearing is introduced. Applications of the developed models to the analysis and optimization of practical cases are discussed. Main outcomes from the work are summarized in the following. 1. Rapid Thermal Predictions: Streamline Upwind Petrov-Galerkin FE Method  The first goal was to formulate a THD model in a robust, efficient and as general as possible form. The streamline upwind Petrov-Galerkin finite element (SUPG-FE) method for three dimensional THD lubrication modeling is developed using Cartesian coordinates. It is shown to be effective in slider as well as journal bearings with some minor modifications. In summary the streamline upwind Petrov-Galerkin FE method can be effectively used to model journal bearings using Cartesian coordinates. The model can be refined, especially the boundary conditions in the groove area, to give more accurate  143  results. The same model can also be used for slider bearings with less implementation requirements. A more detailed model procedure is given in sections 2.4 and 3.3. 2. Slider Bearing Thermal Predictions  Some important observations on slider bearings are noted. An aspect ratio between 2 and 1 is ideal in terms of leakage flow and minimum peak temperatures. A unique way of identifying this region is illustrated where the differential of the flow ∂Q as a function of aspect ratio  ∂Q reduces rapidly to minimum before increasing again and becoming ∂( L / B )  constant at L / B ratios higher than 10. The above findings along with the flow gradient method (FGM), uniquely defined in this work, helps optimize slider bearing tilt ratio for minimum temperature rise. 3. Development of Journal Bearing Template  A template for modeling THD lubrication in journal bearings is shown to be effective in modeling thermal effects, leakage flow, and cavitation effects. The model as presented has no special consideration of thermal boundary conditions at the rupture interface. Although the approach proposed by Boncompain [15] is used in the rupture zone, focus on the dynamics of velocity gradients ∂v/∂y and ∂w/∂y near the rupture and reformation interfaces are critically analyzed and appear to explain the rapid decline in temperature. Thermal predictions are within ten percent of measured values. The model is robust and converges in two or three iterations making it ideal for iterative applications such as in validating optimization equations developed in the second part of the research. 4. Empirically Derived Optimization Model (EDOM)  A novel set of optimization equations derived from first principles is proposed. The set of equations require no complex integrations and can be easily implemented using standard  144  optimization tools such as Matlab. Optimization equations for temperature increase ΔT, power loss PL, and leakage flow QL are developed; constants and other values in the equations are determined empirically for the range of bearings under consideration. The set of equations are well-designed and easy to use while giving accurate results. Comparisons with models in literature show that the equations consistently give better predictions. The approach used by Jang and Khonsari [55] whereby temperature rise parameters are derived and empirically verified is expanded on and used to formulate temperature increase ΔT(X) model. Such an approach is used for the first time in journal bearing optimization. It presents an alternative to Hirani and Suh’s model [38] which requires a finite difference mass conserving algorithm to determine temperature increase and integration techniques to determine leakage flow and power loss. Once empirical values in the temperature increase ΔT(X), leakage QL, and power loss PL are determined, the optimization model can be used for the range of bearings under consideration to find optimum values of the design variables (X). This work presents the developed equations as an optimization model. 5. Optimization using Weight/Scaling Parameters and Pareto Optimal Fronts  The developed set of equations are used to minimize a multiobjective function optimization model using the weight/scaling factors and Pareto optimal fronts and results are compared to Hashimoto’s [34]. It is clear Hashimoto’s model is for short bearings with L/D ratios up to 0.5 and is not accurate for medium range bearings (0.4<L/D<1.0) under consideration here. In the Pareto optimal front approach, posterior articulation of the weights is used whereby the designer initially generates a number of non-inferior (a set of equally efficient) 145  solutions from which a final decision is made on any one solution. Advantages and disadvantages of each method are given, however Pareto charts appear to be preferable as they give a set of non-inferior, non-dominated solutions and final choice depends on the merits of a particular design.  6.2  Recommendations for Future Work  1. Journal Bearing Template  The developed SUPG FE model, although effective in predicting thermal effects, should be extended to thermo-elastohydrodynamic (TEHD) lubrication where expansion of the shaft and bush are taken into account. This will give a better prediction of the load carrying capacity of the bearing. Prediction of the clearance ratio C could be adjusted accordingly. The journal bearing template developed here could be tested on transient phenomenon such as the squeeze effect and thermal predictions during start-up and shut down. This is an important area as it is also linked to the study of bearing time to seizure and how this can be avoided. 2. Journal Bearing System Optimum Design  The model presented here focuses on individual journal bearings. A possible expansion of the work might involve looking at the ‘system’ of a shaft and two journal bearings at the end. Such a task would the overcome the qualifying statement in section 4.3.5 that other factors may influence the choice of design variable set (X) such as the shaft stiffness and misalignment problems that may be encountered.  146  3. Expansion of the Empirically Derived Optimization Model  The optimization model developed here works well for medium range journal bearings. An extension of the model could be developed for short and long range bearings. This would be very useful to design engineers who use bearings in the two ranges and also make the EDOM applicable to the full range of aspect ratios.  147  Bibliography  1. Khonsari, M.M., Wang, S.H., “Note on transient thermo-hydrodynamic lubrication effects in lubricating films”, STLE Tribol. Trans., 35, pp. 177-183, 1992. 2. Pinkus, O., Thermal Aspects of fluid film tribology, ASME Press, New York, 1990. 3. Ettles, C.M.M., Cameron, A., “Considerations of oil flow across a bearing groove”, Trans. ASME J. Lub. Tech., pp. 312-319, 1968.  4. Ettles, C.M.M., “Solutions for flow in a bearing groove”, Proc. Instn. Mech. Engrs., 182, p.120, 1967-8.  5. Kosasih, P.B., Tieu, A.K., “An investigation into the thermal mixing in journal bearings”, Proc. Inst. Mech. Engrs. Part J: J. Eng. Tribol., 218, pp. 379-389, 2004. 6. Costa, I., Miranda, A.S., Fillon, M., Claro, J.C.P., “An analysis of the influence of oil supply conditions on the thermo-hydrodynamic performance of a single-groove journal bearing”, Proc. Inst. Mech. Engrs. Part J: J. Eng. Tribol., 217 (J2), pp. 133144, 2003. 7. Rodkiewicz, C.M., Sinha, P., “On the lubrication theory: a mechanism responsible for generation of the parallel bearing load capacity”, Trans. ASME J. Tribol., 105, pp. 584-590, 1993. 8. Schumack, M.R., “Application of the pseudo-spectral method to thermohydrodynamic lubrication”, Int. J. Numer. Methods. Fluids, 23, pp. 1145-1161, 1996. 9. Ratish Kumar, B.V., Srinivasa Rao, P., Sinha, P., “A SUPG/FEM analysis of thermal effects on load carrying capacity in slider bearing lubrication”, Numer. Heat Transfer Part A, 38(3), pp. 305-328, 2000.  148  10. Ratish Kumar, B.V., Srinivasa Rao, P., Sinha, P., “A numerical study of performance of a slider bearing with heat conduction to pad”, Fin. Elem. Analysis Design, 37, pp. 533-547, 2001. 11. Brooks, A.N., Hughes, T.J.R., “Streamline upwind/Petrov-Garlerkin formulations for convective dominated flows with particular emphasis on the incompressible NavierStokes equations”, Comp. Methods Appl. Mech. Eng., 32, pp. 199-259, 1982. 12. Grygiel, J-M., Tanguy, P.A., “Finite element solution for advection-dominated thermal flows”, Comput. Methods Appl. Mech. Eng., 93(3), pp. 277-289, 1991. 13. Zengeya, M., Gadala, M., “Viscous dissipation effects in hydrodynamic lubrication”, Proc. Zim. Inst. Engrs., 4(1) , pp. 5-16, 2006.  14. Ezzat, H.A., Rohde, S.M., “A study of the thermo-hydrodynamic performance of finite slider bearings”, Trans. ASME J. Lub. Tech., 95(3), pp 298-307, 1973. 15. Boncompain, R., Fillon, M., Frene, J., “Analysis of thermal effects in hydrodynamic bearings”, Trans. ASME J. Tribol., 108, pp. 219-224, 1986. 16. Gero, L.R., Ettles, C.M., “A three-dimensional thermo-hydrodynamic finite element scheme for fluid film bearings”, ASLE Tribol. Trans., 31(2), pp 182-191, 1988. 17. Hahn, E.J., Kettleborough, C.F., “Solution for the pressure and temperature in an infinite slider bearing of arbitrary profile”, Trans. ASME J. Lub. Tech., 89, p 445, 1967. 18. Kucinschi, B.R., Dewitt, K.J., Pascovici, M.D., “Thermoelastohydrodynamic (TEHD) analysis of a grooved thrust bearing”, Trans. ASME J. Tribol., 126(2), pp. 267-274, 2004.  149  19. Paranjpe, R. S., Goenka, P.K., “Analysis of crankshaft bearings using a mass conserving algorithm”, STLE Tribol. Trans., 33(3), pp 333-344, 1990. 20. Paranjpe, R.S., Han, T., “A finite volume analysis of the thermo-hydrodynamic performance of finite journal bearings”, Trans. ASME J. Tribol., 112, pp 557-566, 1990. 21. Paranjpe, R.S., Han, T., “A study of the thermo-hydrodynamic performance of steadily loaded journal bearings”, STLE Tribol. Trans., 37(4), pp 679-690, 1994. 22. Stachowiak, G.W., Batchelor, A.W., Engineering Tribology, 2nd Ed., Butterworth, p. 216, 2001. 23. Booker, J.F., Huebner, K.H., “Application of the finite element methods to lubrication: an engineering approach”, Trans. ASME J. Lub. Tech., pp 313-323, 1972. 24. Goenka, P.K., “Dynamically loaded journal bearings: finite element method analysis”, Trans. ASME J. Tribol., 106, pp 429-439, 1984. 25. Labouff, G.A., Booker, J.F., “Dynamically loaded journal bearings: a finite element treatment for rigid and elastic surfaces”, Trans. ASME J. Tribol., 107(4), pp 505-515, 1985. 26. Bayada, G., Chambat, M., El Alaoui, “Variational formulations and finite element algorithms for cavitation problems”, Trans. ASME J. Tribol., 112, pp 398-403, 1990. 27. Kucinschi,  B-R.,  Fillon,  M.,  Frene,  J.,  Pascovici,  M.D.,  “A  transient  thermoelastohydrodynamic study of steadily loaded plain journal bearings using finite element method analysis”, Trans. ASME J. Tribol., 122(1), pp 219-226, 2000. 28. Elrod, H.G., “A cavitation algorithm”, Trans. ASME J. Lub. Tech., 103(3), pp. 350354, 1981.  150  29. Kumar, A., Booker, J.F., “A finite element cavitation algorithm”, Trans. ASME J. Tribol., 113, pp 276-286, 1991.  30. Kumar, A., Booker, J.F., “A mass and energy conserving finite element lubrication algorithm”, Trans. ASME J. Tribol., 116, pp 667-671, 1994. 31. Optasanu, V., Bonneau, D., “Finite element mass-conserving cavitation algorithm in pure squeeze motion”, Trans. ASME J. Tribol., 122, pp 162-169, 2000. 32. Colynuck, A.J., Medley, J.B., “Comparison of two finite difference methods for the numerical analysis of thermo-hydrodynamic lubrication”, STLE Tribol. Trans., 32, pp. 346-356, 1989. 33. Shigley, J.E., Mischke, C.R., Budynas, R.G., Mechanical Engineering Design, 7th Ed., McGraw-Hill, New York, p. 619, 2004. 34. Hashimoto H., “Optimum design of high-speed, short journal bearings by mathematical programming”, STLE Tribol. Trans., 40(2), pp. 283-293, 1997. 35. Hashimoto H., Matsumoto K., “Improvement of operating characteristics of highspeed hydrodynamic journal bearings by optimum design: part I – Formulation of methodology and its application to elliptical bearing design”,  Trans. ASME J.  Tribol., 123, pp. 305-312, 2001.  36. Song J-D, Yang B-S, Choi B-G, Kim H-J, “Optimum design of short journal bearings by enhanced life optimization algorithm”, Tribol. Int., 38, pp. 403-412, 2005. 37. Yang B-S., Lee Y-H., Choi B-K., Kim H-J., “Optimum design of short journal bearings by artificial life algorithm”, Tribol. Int., 34, pp. 427-435, 2001. 38. Hirani H, Suh NP, “Journal bearing design using multiobjective genetic algorithm and axiomatic design approaches”, Tribol. Int., 38, pp. 481-491, 2005.  151  39. Suh NP, “Axiomatic design: advances and applications”, Oxford University Press, 2001. 40. Hirani, H., “Multiobjective optimization of journal bearing using the pareto optimality concept”, Proc. Inst. Mech. Engrs. Part J: J. Eng. Tribol., 218(4), pp. 323-336. 41. Hirani, H., “Multiobjective optimization of journal bearing using mass conserving and genetic algorithms”, Proc. Inst. Mech. Engrs. Part J: J. Eng. Tribol., 219(3), pp. 235-248, 2005. 42. Deb K, Multiobjective optimization using evolutionary algorithms, John Wiley & Sons, New York, 2001. 43. Hashimoto, H., “Optimization of oil flow rate and oil film temperature rise in high speed hydrodynamic journal bearing”, in: Dowson D., editor, Tribology for energy conservation, Amsterdam, Elsevier, pp. 205-210, 1998.  44. Powell, M.J.D., “A fast algorithm for non-linearly constrained optimization calculations”, GA Watson Ed., Lecture Notes in Mathematics, Springer Verlag, 630, 1978. 45. Goldberg, D.E., “Genetic algorithm in search, optimization and machine learning”, Addison Wesley Publishing Company, 1989. 46. Holland, J.H., “Adaptation in natural and artificial systems”, University of Michigan Press, Michigan, 1975. 47. DeJong, K., “The analysis and behavior of genetic adaptive systems”, PhD Dissertation, University of Michigan, Michigan, 1975.  152  48. Homaifar, A., Qi, C.X., Lai, S.H., “Constrained optimization via genetic algorithms, simulations”, Engineering Optimization, 62, pp. 242-254, 1994. 49. Louis, S.J., Zhoo, F., and Zeng, X., “Flaw detection and configuration with genetic algorithms”: in Evolutionary Algorithms in engineering applications, SpringerVerlag, 1997. 50. Saruhan, H., “Design of a three-lobe journal bearing using a genetic algorithm”, J. Naval Sci. Eng., 2(2), pp. 87-103, 2004.  51. Genetic Algorithm and Direct Search Toolbox User’s Guide for use with Matlab, The Mathworks Inc. 2004. 52. Zengeya, M., Gadala, M., Segal, G., “Three-dimensional thermal field in slider bearings”, Proc. World Tribology Congress III, Washington DC, p. 30, September 12-16 2005. 53. Zengeya, M., Gadala, M., Segal, G., “Three-Dimensional Modeling of Thermohydrodynamic Lubrication in Slider Bearings Using Streamline Upwind PetrovGalerkin Method”, Num. Heat Transfer, Part A, 49(10), pp. 947-968, 2006. 54. Dowson, D., “A generalized Reynolds equation for fluid-film lubrication”, Int. J. Mech. Sci., 4, pp 159-170, 1962.  55. Jang, J.Y., Khonsari, M.M., “Design of bearings on the basis of thermodrodynamic analysis”, Proc. IMechE J. Eng. Trib., 218, pp 355-363, 2004. 56. Segal, G. (2003), Sepran Standard Problems, Ingenieursbureau SEPRA, Den Haag, The Netherlands. 57. Ezzat, H.A., Rohde, S.M., “A study of the thermo-hydrodynamic performance of finite slider bearings”, Trans. ASME J. Lub. Tech., 95(3), pp. 298-307, 1973.  153  58. Boncompain, R., Frene, J., “Thermo-hydrodynamic analysis of finite journal bearings – static and dynamic characteristics”: in Proc. 6th Leeds-Lyon Symposium on Tribology: Thermal Effects in Tribology, pp. 33-41, 1980.  59. Dowson, D., Hudson, J., Hunter, B., March, C.N., “An experimental investigation of the thermal equilibrium of steadily loaded journal bearings”, Proc. IMechE. Eng., 181 (3B), pp 70-80, 1966. 60. Reddi, M.M., “Finite element solution of the incompressible lubrication problem”, Trans. ASME J. Lub. Tech., 69, pp 524-531, 1969.  61. Hahn, E.J., Kettleborough, C.F., “Solution for the pressure and temperature in an infinite slider bearing of arbitrary profile”, Trans. ASME J. Lub. Tech., 89, p. 445, 1967. 62. Ingenieursbureau SEPRA, Park Nabij 3, 2491 EG Den Haag, The Netherlands. 63. Zengeya, M., Gadala, M., Segal, G. (2007). Hydrodynamic and thermal behavior of journal bearings using streamline upwind Petrov-Galerkin FEM, STLE Tribology Trans., 50(2), 227-247.  64. Hamrock, B.J., “Fundamentals of fluid film lubrication”, 2nd Edition, Marcel Dekker, New York, p168, 2004. 65. Lebeck, A.O., “Parallel sliding load support in the mixed friction regime”, Trans. ASME J. Trib., 109, pp. 189-195, 1987b.  66. Arnell, R.D., Davies, P.B., Halling, J., and Whomes, T.L., Tribology, principles and design applications, Springer-Verlag, New York, 1991.  67. Ezzat, H.A., Rohde, S.M., “A study of the thermo-hydrodynamic performance of finite slider bearings”, Trans. ASME J. Lub. Tech., 95(3), pp 298-307, 1973.  154  68. Khonsari, M.M., Booser, E.R., Applied tribology: bearing design and lubrication, John Wiley, New York, p 215, 2001. 69. Costa, L., Fillon, M., Miranda, A.S., Claro, J.C.P., “An experimental investigation of the effect of groove location and supply pressure on the THD performance of a steadily loaded journal bearing”, ASME J. Trib., 122, pp 127-132, 2000. 70. Martin, F.A., “Oil flow in plain steadily loaded journal bearings: realistic predictions using rapid techniques”, Proc. IMechE J. Eng. Tribol., 212, 1998, 413-425. 71. Kelly, D.W., Nakazawa, S., Zienckiewicz, O.C., “A note on upwinding and anisotropic balancing dissipation in finite element approximations to convective diffusion problems”, Int. J. Num. Meth. Eng., 15, pp. 1705-1711, 1980. 72. J. Mitsui, “A study of thermo-hydrodynamic lubrication in a circular journal bearing”, Tribol. Int., 20(6), 331-341, 1987.  73. Pierre, I., Fillon, M., “Influence of geometrical parameters and operating conditions on the thermo-hydrodynamic behavior of plain journal bearings”, Proc. IMechE J. Eng. Tribol., 214, 2000, 445-457.  74. Khonsari, M.M., Jang, J.Y., Fillon, M., “On the generalization of thermohydrodynamic analysis for journal bearings”, Trans. ASME J. Trib., 118, pp. 571-579, 1996. 75. Ferron, J., Frene, J., Boncompain, R., “A study of the thermo-hydrodynamic performance of a plain journal bearing: comparison between theory and experiments”, Trans. ASME J. Lub. Tech., 105, 422-428, 1983. 76. Zengeya, M, Gadala, M., “Optimization of journal bearings using a hybrid scheme”, IMechE Part J: J. Eng. Tribol., 221(5), 591-607, 2007.  155  77. Khonsari M. M., Booser E.R., "Parasitic Power Losses in Hydrodynamic Bearings". Machinery Lubrication, March 2006.  78. Optimization Toolbox for use with Matlab, User’s Guide Version 2, The Mathworks Inc., 2003. 79. Claro, J.C.P., and Miranda, A.A.S., “Analysis of hydrodynamic journal bearings considering lubricant supply conditions”, IMechE Part C: J. Mech. Eng. Sci., 207, 93-101, 1993. 80. Dowson, D., Taylor, C.M., Miranda, A.A.S., “The prediction of liquid film journal bearing performance with a consideration of lubricant film reformation Part 1: theoretical results”, Proc. IMechE J. Eng. Trib., 199(C2), 95-102, 1985. 81. Dowson, D., Taylor, C.M., Miranda, A.A.S. “The prediction of liquid film journal bearing performance with a consideration of lubricant film reformation Part 2: experimental results”, Proc. IMechE J. Eng. Trib., 199(C2), 103-111, 1985.  156  APPENDICES APPENDIX A A.1 Velocity fields u, v, and w  The equilibrium condition for forces acting in the x-direction gives ∂p ∂ ⎛ ∂u ⎞ ⎟=0 − ⎜μ ∂x ∂y ⎜⎝ ∂y ⎟⎠  (A-1)  where u is the velocity of the lubricant in the direction of sliding. The boundary conditions are u=0 at y=h and u=U at y=0. For the one-dimension bearing considered here the value of the lubricant flow is constant along the bearing assuming no side leakage; this fact is used to calculate v. Algebraically this means that the lubricant  ( ) h  velocity across the film thickness is constant, that is ∂ ∂x ∫ udy = 0 . Equation A-1 can be 0  integrated in the y-direction to get v, the component of u in that direction . u( y) =  dp dx  ∫  y  0  y  μ  dy + A  ∫  y  0  dy  μ  +U2  (A-2)  where A is given by ⎛ ⎜⎜ − dp dx A= ⎝  ⎞  h  ∫ ( y μ )dy ⎟⎟⎠ − U ∫ dy μ 0  2  (A-3)  h  0  and U the slider velocity in m/s. Non-dimensional value becomes  /  . Similar  equilibrium condition for forces acting in the y-direction is  ∂p ∂ ⎛ ∂w ⎞ ⎟=0 − ⎜μ ∂z ∂y ⎜⎝ ∂y ⎟⎠  (A-4)  with boundary conditions w=0 at y=0 and w=0 at y=h. Equation A-4 is integrated with respect to y to give  157  ⎡ dp ⎢ y ( y μ )dy − w= dz ⎢ ∫0 ⎣⎢  (∫ dy μ )(∫ ( y μ )dy )⎤⎥ y  h  0  0  h  ∫ dy 0  μ  ⎥ ⎦⎥  (A-5)  Since pressure is constant through the film thickness, ∂p ∂y = 0 and the continuity of flow requires that ∂u ∂v ∂w + + =0 ∂x ∂y ∂z  (A-6)  The continuity equation is more conveniently solved by differentiating with respect to y giving ∂ 2 v ∂ ⎛ ∂u ∂w ⎞ + + ⎜ ⎟=0 ∂y 2 ∂ y ⎝ ∂x ∂z ⎠  with boundary conditions u=w=0 at y=h and u=U, w=0 at y=0  158  (A-7)  APPENDIX B  The term (1+εcosϕ)3 in equation 5.12, referred to as (hg/C)3 in Martin’s model, relates to the case where the oil groove is line with the load. The same term applies for an inlet hole of diameter DH in line with the load but fg in this case becomes ⎛D ⎞ f g = 0.675⎜ H + 0.4 ⎟ ⎝ L ⎠  1.75  0.1<DH<0.25  (B.1)  and '  Q L ,tot = QmS Q 1p− S  '  (B.2)  where S’=0.6. Another important case is when the oil groove is at the maximum film thickness position hmax where (hg/C)3 becomes  ⎛ hg ⎜⎜ ⎝C  3  ⎞ ⎟⎟ = (1 + ε )3 ⎠  (B.3)  equation B1 for fg is valid, and total leakage becomes  QL ,tot  ⎛ Lg = Qm ⎜⎜ ⎝ L  ⎞ ⎟⎟ ⎠  m  ⎛Q ⎞ where m = 0.27⎜ L Q ⎟ p ⎠ ⎝  0.3<Lg/L<0.8  (B.4)  0.27  . The proposed model here uses equation 10 to determine the  flow due to film pressure QL while Qp is evaluated using equation 13 with the appropriate value of (hg/C)3. The power value to eccentricity ratio in equation B4 applies when the hydrodynamic film starts at hmax and L/D<1.0, 0.3≤λ≤0.9.  159  APPENDIX C Optimization Procedure  The procedure presented here will help designers implement the optimization model proposed. This can be easily done in Matlab.  A1. Sequential Quadratic Programming (SQP)  This requires three main components namely the main program, objective function, and non-linear constraint m-files. The call to the optimization subroutine ‘fmincon’ is done in the main program while computations in part B are done in the objective function m-file and part C computes the non-linear constraints. Define global variables that will be passed between the three m-files (ΔT, ε, Ns). The three m-files are attached at the end of the appendix. While care was taken in preparing these files, no guarantee is given that the codes are error free. A brief description of the main parts follows. Create a folder and copy the three files. Run the main program which will call the other m-files.  Part A: Main Program -Define upper and lower limits for design variable X T = (C , λ , μ , p s ) and provide these as lb and ub in the call to ‘fmincon’ -Provide program with initial guess for design variable x0 -Set options and call ‘fmincon’ (this calls both Part B and C)  160  Part B: Objective Function -Set all the problem parameters here (since the hybrid technique will use same objective function) For Ns from Ns1 to Ns max Step 1: Compute journal surface velocity U 2: Use μI to compute Sommerfeld number S 3: Find eccentricity ratio ε 4: Compute the temperatures Tmax and Tshaft and temperature increase ΔT 5: Update the viscosity μ 6: Go to step 2 and re-compute S until convergence. 7: Compute the side leakage QL 8: Compute the objective function F(X) given the constraints gi(X) 9: Record F(X) 10. Loop and update Ns  Part C: Non-linear Constraint Computation -Define minimum allowable hmin, ΔTa, ωcr ,and pa -Compute pmax, ω, and ωcr -Define non-linear constraints g1, g2, g3, and g4  C.2 Matlab M-files  Main Program function [X,FVAL,REASON,OUTPUT,POPULATION,SCORES] = dowprob1 %% This is an auto generated M file to do optimization with the Genetic Algorithm and % Direct Search Toolbox. Use GAOPTIMSET for default GA options structure. %%Fitness function fitnessFunction = @dowtaylorobj; %%Number of Variables  161  nvars = 5 ; %Linear inequality constraints Aineq = []; Bineq = []; %Linear equality constraints Aeq = []; Beq = []; %Bounds LB = [4e-005 0.4 0.01 50000 1.5708 ]; UB = [0.0003 1 0.03 350000 1.5708 ]; %Nonlinear constraints nonlconFunction = []; %Start with default options options = gaoptimset; %%Modify some parameters options = gaoptimset(options,'PopInitRange' ,[-2.5 ; 2.5 ]); options = gaoptimset(options,'PopulationSize' ,100); options = gaoptimset(options,'MutationFcn' ,{ @mutationgaussian 1 1 }); options = gaoptimset(options,'Display' ,'off'); options = gaoptimset(options,'PlotFcns' ,{ @gaplotbestf @gaplotbestindiv }); options = gaoptimset(options,'HybridFcn' ,{ @fminsearch [] }); %%Run GA [X,FVAL,REASON,OUTPUT,POPULATION,SCORES] = ga(fitnessFunction,nvars,Aineq,Bineq,Aeq,Beq,LB,UB,nonlconFunction,opti ons);  Objective Function function [ret,f1,f2] = objective(x) global ns rho c_p D W e_0 deltaT Q_tot rfric1 p_max p_s f1 f2 a1 b1 c1 omega1 omega_cr tht_mn; global tht_mx p_smin p_smax; n_s = 3; n_smx =50; rho = 860; c_p = 2.0e3; mu_i = 0.0358; T_i1 = 313; T_i = T_i1 - 273; k_f = 0.13;  % Journal speed, 40-240 rps  alpha = 0.756e-7; beta = 0.0414;  % Thermal diffusivity of lubricant, m^2/s % Temperature-viscosity coefficient, K^(-  % % % % % %  Density of fluid, kg/m3 Specific heat of lubricant Inlet viscosity of oil, Pa.s Inlet oil temperature, K Inlet temp in C Thermal conductivity of lubricant, W/(m  K)  1) L1 = 0.08; L_g1 = 0.07; L_Lg1 = L1/L_g1; w_g = 12e-3; D = 0.1; W = 20e3; C_mn = 40.0e-6;  % Bearing length % Length of groove % Ratio of groove/bearing length % Width of groove % Journal diameter % Load on bearing, N % Radial clearance  162  C_mx = 300.0e-6; lambda_mn = 0.20; lambda_mx = 1.0; mu_mn = 1.0e-3; mu_mx = 4.0e-2; p_smin = 50e3; p_smax = 350e3; beta1 = 1; beta2 = 1e5; alpha1 = 1; alpha2 =alpha1;  % L/D ratio % Viscosity % Supply pressure % Scaling factor % Weighting factor  % Nonlinear objective function %x0 = [50e-6 0.6 2e-2];  %  for ns = n_s:n_smx, % Define the variables C=x(1); La=x(2); mu=x(3); p_s=x(4); theta = x(5); L = La*D; L_g = L_g1; R_C = D/(2*C); ps_bar = p_s*(1/R_C)^2/(mu*2*pi*ns); U = ns*2*pi*D/2; Re = rho*C*U/mu; S = ns*mu*D*L*(D/(2*C))^2/W; e_0 = -0.36667+0.49737/sqrt(La)-0.20986*log(S); if (e_0 > 1), e_0 = 0.6121; return; end;  % Compute side leakage Q_Leak Q_leak1 = exp(0.20-0.5*La^1.5+1.05*log(e_0)); Q_bar4 = (1+1.6*ps_bar)*(1+0.17*theta)*(1+0.3*w_g/D)*(1.15*L_g/L)*Q_leak1; Q_leak4 = Q_bar4*pi*ns*D*C*L; Q_leak = Q_leak1*pi*ns*D*C*L; % find leakage due to supply pressure z = 4*(1 + La*(1 - 1.25*e_0^0.18)); attit = atan((pi/(z*e_0))*sqrt(1-e_0^2)); f_g = (1.25-0.25*(L_g/L))/(3*(L/L_g-1)^(0.33333))+... (w_g/D)/(3*La*(1-L_g/L)); h_g = (1 + e_0*cos(attit))^3; Q_p =f_g*h_g*p_s*C^3/mu_i; % Leakage due to pressure Qm = Q_leak + Q_p - 0.3*sqrt(Q_leak*Q_p); % Datum flow rate S_prime = 0.7*(L_g/L)^0.7 + 0.4; Q_tot = Qm^S_prime*Q_p^(1-S_prime); % Total leakage  % % %  %Compute temperature increase kappa1 = alpha*mu_i*beta*U*(D/(2*C))^2/(k_f*D/2); kappa2 = U*sqrt(mu_i*beta/k_f); % Determine parameters a, b, and c if (e_0 > 0.1) && (e_0 < 0.3), a1 = 1.58105; b1 = -0.14250; c1 = -1.04050; end; if (e_0 >= 0.1) && (e_0 <= 0.7),  163  a1 = 1.5917; b1 = -0.13962; end; if (e_0 > 0.7) && (e_0 <= 0.9), a1 = 1.63973; b1 = -0.10627; end;  c1 = -1.24999;  c1 = -1.45216;  % Compute T_max T_max1 = exp(a1 + b1/kappa1^0.5 + c1/kappa2^0.5); T_max = (T_max1/beta + T_i) + 273; deltaT = T_max - T_i1; % Compute Power loss rfrict = exp(1.9994+2.73622/sqrt(La)-4.97838*sqrt(e_0)); rfric1 = (10.25*theta^2+0.23*abs(theta))*(1+0.6*ps_bar)*(rfrict); torqueM = rfrict/R_C*W*D/2; torq2 = rfric1/R_C*W*D/2; P_Loss = torqueM*2*pi*ns; Ploss2 = torq2*2*pi*ns; % Recompute mu x(3) = mu_i*exp(-beta*deltaT); % Objective function f1 = alpha1*beta1*Ploss2; function f2 = alpha2*beta2*Q_leak4; ret = f1 + f2;  % Friction parameter obj  end  Non-linear Constraints M-file function [c, ceq] = nonlin_cons(x) global ns rho c_p D W e_0 deltaT Q_tot rfric1 p_max p_s f1 f2 a1 b1 c1 omega1 omega_cr tht_mn; global tht_mx p_smin p_smax; h_a = 5e-6; p_a = 30e6; deltaT_a = 70;  % Allowable minimum film thickness % Allowable maximum pressure, Pa % Allowable temperature rise, K  p_max1 = exp(2.25062-0.39491/x(2)+4.77721*e_0^3); p_max = p_max1*x(3)*(D/2)^2*ns/x(1)^2 + p_s; omega1 = 2*pi*ns; omega_cr = sqrt(9.81/x(1))*(0.0584*exp(6.99*e_0^2.07)-1.318*e_0+2.873); g7 = h_a -x(1)*(1-e_0); g8 = deltaT - deltaT_a; %g9 = omega1 - omega_cr; g10 = p_max - p_a; c = [g7 g8 g10]; % nonlinear equality constraints returned as vector ceq ceq = []; % there are no nonlinear equality % constraints in  164  C3. Hybrid Scheme  The genetic algorithm toolbox with the Hybrid option was used. The fitness function is the objective function defined in section C1 and non-linear constraints are the same as in C1. Other settings are as shown in Table D5.  165  APPENDIX D Tables  Table D.1 Input data from Dowson [59] and Song [36] bearings (empty cells mean unspecified values). Input Parameter  Dowson  Song  Bearing Diameter D, m Bearing width W, m Bush outer diameter, m  0.102 0.07625 0.2 8.01 m/s 11000 0.76 800 4.76 x 10-3 m 0.067 m 2000, 380, 490 0.13, 50, 65 W/m K 100  0.1 10000 4190 -  200  -  308 K 867 kg/m3 0.03 Pa.s 0.0414 /K  -  Density-temperature coeff.  0.0012 /K  -  Thermal diffusivity (αt=ρcp/kf) Supply pressure ps (gauge)  0.756 x 10-7 m2/s  -  0.276 x 106 Pa  -  Journal Speed (Ns=25rev/s), U  Load, W L/D Radial clearance ratio R/C Groove width Groove length Specific heat cf, cb, cs Conductivities kf, kb, ks Convective heat transfer coefficient H1, W/m2.K Film/bush heat transfer coefficient H2, W/m2.K Inlet temperature Ti Inlet lubricant density ρi Inlet viscosity μi Viscosity-temperature coeff.  β δ  D1.1 Location of thermocouples  Figure D1 shows the location of thermocouples (TC) on the bush and shaft in Dowson’s experimental setup [59]. The temperature of the lubricant at inlet was recorded by two iron-constantan thermocouples fitted in a radial inlet hole located within the bush. One of the thermocouples was located near the outer radius of the bush while the other was  166  within 1.59 mm of the entry to the oil groove. Four outlet temperature thermocouples were fitted in special oil scoops located at the ends of the bush to collect lubricant flow leaving the bearing. For bush temperature measurements, the thermocouples were mounted in groups of four on bronze plugs which fitted closely into 34 radial holes in the bush. The hot junction of each TC was located on the centerline of each plug and spaced at distances of 1.59, 17.5, 33.3, and 49.2 mm from the dead end of the radial hole adjacent to the oil film. The bronze plugs were located in four planes at the axial locations shown in figure D1a. Shaft temperatures were measured using twelve TC’s mounted in plugs similar to the bush setup. The TC plugs were tightly fitted into radial holes in the shaft with their outer ends flush with the shaft surface. The TC’s were located in four planes disposed about the edge of the bearing as shown in figure D1b to enable temperature gradients to be determined. Wires from the thermocouples ran along a hole drilled in the shaft to a slipring assembly that was maintained at a constant temperature by air cooling at the end of the shaft. The shaft surface temperature was recorded by means of a specially constructed thin film platinum resistance thermometer mounted on a glass rod located in a radial hole in the shaft at a distance of 16 mm from the central plane of the bearing. The exposed surface of the glass plug was contoured to conform with the shaft surface. The leads fro the platinum film were carried down the internal bore of the shaft to the slip ring assembly.  167  Fiigure D1. Th hermocouplee location forr pressure annd temperatuure measurem ments [59]  168  Table D.2 Input values for the bearing tested by Boncompain et al. [15] Input Parameter  Value  Bearing Radius R, m Bearing width W, m Bush outer radius R2, m Journal Speed (Ns=33.3rev/s), U Eccentricity, ε Radial clearance at 20°C, C, m-3 Groove angle Groove length Groove Location Specific heats cf, cb, cs Conductivities kf, kb, ks, kg, W/m K Convective heat transfer coefficient , bush/ambient H1, W/m2.K Convective heat transfer coefficient , bush/oil and shaft/oil interfaces H2, W/m2.K Inlet temperature Ti, K Ambient Temperature To, K Inlet lubricant density ρi Inlet viscosity μi Viscosity-temperature coeff. β Density-temperature coeff. δ Supply pressure ps (gauge)  0.050 0.08 0.100 10.5 m/s 0.8 145µm 18° 0.070 m 0° 2000, 380, 490 0.13, 250, 50, 0.025  169  100  500 313 318 860 kg/m3 0.028 Pa.s 0.0414 /K 0.0012 /K 0.070 x 106 Pa  Table D.3 Input values for Costa’s bearing [69] Input Parameter  Value  Bearing Diameter D, m Bearing width W, m Bush outer diameter, m Journal Speed (Ns=50rev/s), U Load, W L/D Radial clearance C, m-6 Groove width Groove length Groove Location Specific heats cf, cb, cs Conductivities kf, kb, ks Convective heat transfer coefficient H, W/m2.K Inlet temperature Ti Inlet lubricant density ρa Inlet viscosity μi Viscosity-temperature coeff. β Density-temperature coeff. δ Supply pressure ps (gauge)  0.100 0.08 0.2 15.7 m/s 8000 N 0.75 93.5 16 x 10-3 m 0.070 m 0° 2000, 380, 490 0.13, 50, 65 W/m K 100  170  308 K 867 kg/m3 0.0358 Pa.s 0.0414 /K 0.0012 /K 0.300 x 106 Pa  Table D.4 Input values for Ferron’s [75] bearing. Input Parameter  Ferron  Bearing Radius R, m Bearing width L, m Bush outer radius R2(assumed), m Radial clearance at 20°C, C, m-6 Groove angle, ° Groove length, m-3 Groove Location, ° Specific heats cf, cb, cs Conductivities kf, kb, ks, kg, W/m K Convective heat transfer coefficient , bush/ambient H1, W/m2.K Convective heat transfer coefficient , bush/oil and shaft/oil interfaces H2, W/m2.K Inlet temperature Ti, K Ambient Temperature To, K Inlet lubricant density ρI, kg/m3 Inlet viscosity μi, Pa.s Viscosity-temperature coeff. β, /°K Supply pressure ps (gauge), kPa  0.05 0.08 0.1 145 70 (assumed) 0 2000, 380, 490 0.13, 250, 50, 0.025  171  80  250 313 313 860 0.0277 0.0414 70  Table D.5 Additional input data to run optimization simulations. Parameter Minimum radial clearance Maximum radial clearance Minimum length to diameter ratio λ Maximum length to diameter ratio λ Minimum lubricant viscosity, Pa.s Maximum lubricant viscosity, Pa.s Minimum allowable film thickness Maximum allowable film pressure Allowable temperature rise Scaling (normalizing) factor Weighting factor Load, kN Population Size (Hybrid Technique) Population Type Initial Search Range Fitness Scaling Selection Function Reproduction: Cross-over Fraction Mutation Function Cross-over Function Migration Direction  Dowson Cmin=40μm  Song Cmin=40μm  Cmax=300μm  Cmax=300μm  λmin=0.4  λmin=0.2  λmax=1.0  λmax=0.6  μmin=0.01  μmin=0.001  μmax=0.03  μmax=0.01  hmin=10 μm  hmin=10 μm  pmax=10 MPa  pmax=10 MPa  ΔTa=70 °C  ΔTa=70 °C  β1=1, β2=105  β1=1, β2=105  α1=1, α2=α1/5  α1=1, α2=α1/5  10, 20 100  10, 20 20  Double Vector -2.5 to 2.5 Rank Stochastic Uniform 0.8 Gaussian Scattered Forward  172  Table D6 Input data for Dowson [80] and Claro [79] bearings  Dowson et al.  Claro  0.0635 0.0635  0.05 0.05  138.8 0.026 12x10-3 0 2000, (840), (490) 0.13, (1.05), 50, 0.025 313 (298)  125 0.04 10x10-3 +90° (2000, 380, 490) 0.13, 250, 50, 0.025 313 (298)  867.5 0.0224  860 0.075  (0.0414)  (0.0414)  0,170 700  8-60 824  Input Parameter  Bearing diameter D, m Bearing width L, m Radial clearance at 20°C, C, μm Groove length Lg, m Groove width wg, m Groove Location, ° Specific heats cf, cb, cs Conductivities kf, kb, ks, kg, W/m K Inlet temperature Ti, K Ambient Temperature To, K Inlet lubricant density ρI, kg/m3 Inlet viscosity μi, Pa.s Viscosity-temperature coeff. β, /°K Supply pressure ps (gauge), kPa Load W, N  173  APPENDIX E Verification Examples E1  Introduction  This section presents some further validation examples of the streamline upwind PetrovGalerkin finite element method developed. The classical method is compared to the Galerkin and SUPG methods, namely the classical upwind scheme defined in section 2.5, the doubly asymptotic, and discontinuity capturing methods. Brief definitions of each follow. E1.1  Upwinding Methods  If in the energy equation 2.11 the convective part  u·  dominates the diffusive part  , an improvement of the accuracy may be possible by applying the upwinding method. It should be noted that upwinding does not always improve the accuracy and moreover is more computationally expensive in building the matrices. The SUPG method used in this work is based on the method proposed by Brooks and Hughes [11]. Essentially the method introduces an extra term to the standard Galerkin equation  Ω  (E.1)  where e is the element, DT is the differential equation (energy equation in this case), f the source term, and  the weighting function that contains two parts, the Galerkin weighting  function and ζ defined in equation 2.30. The parameter defining the type of upwinding ξ in equation 2.30 defines the type of upwinding. For the definition of ζ we introduce  η=uTΛu where Λ is the matrix with elements Λij and ν becomes  174  u ∆  η  . The type of upwinding  1. Classical SUPG upwinding: ζ=1 which is used in all the results in thesis. ν  2. Doubly asymptotic approximation: ξ  3≤ν≤3           ν           |ν|>3  3. Discontinuity capturing type 1 from Hughes et al. [11] where a pseudo velocity u  u· c  where c is the degree of freedom (solution to be computed) is defined  and the upwind basis function ζ=(τ1u+τ2u )· ∆  the basis functions and  and  where ∇Ni is the divergence of ∆  . For this case ξ is chosen  according to the doubly asymptotic approximation.  (E.2)  More upwinding schemes appear in the work of Brooks and Hughes [11], however only these three will be analyzed and results compared to the Galerkin method. The first example analyses a convection-diffusion problem while the second compares SUPG and Galerkin on lubrication data.  E.1.2 Convection Skew to the Mesh E.1.2.1 Problem Statement  Consider a convection-diffusion problem defined by the general equation u·  (E.3)  where D is the diffusivity and c is the concentration of a molecular species. We assume steady state conditions (∂c/∂t=0) and no source term (fs=0). We consider the flow in the square shown in figure E1 with diffusivity D=10-6 with unidirectional and constant flow u  1 (convection dominates the flow) .  Boundary conditions: Lower boundary (y=0):  Dirichlet (c=1)  175  Dirichlet (c=1 for y≤0.2 and c=0 for y>0.2)  Left boundary (x=0): Natural boundary condition,  0 imposed on all other boundaries  Flow Conditions: The flow is unidirectional and constant ( u  1)  The angle of the flow α is an input variable chosen to be 45°. y  1  c=0  α  0.2 c=1 0  c=1  x 1  0  Figure E1. Region for skew convection  Exact Solution: The exact solution is equal to one in the region starting with lower boundary condition c=1 and following a straight line with the angle of flow. Obviously the value of c lies between o and 1 (0≤c≤1)  E.1.2.2 Solution Procedure  The diffusivity D is much less than flow velocity u thus the problem is convectiondominated. The behavior of the standard Galerkin and the three upwind methods is compared. Linear triangular elements are used. Figure E2 shows the mesh.  176  Figure E2. Mesh for convection-diffusion problem  E.1.2.4 Simulation Results  Appearance of wiggles in the solution is evident in the minimum and maximum values. Table E1 compares results from the four methods.  Table E1. Minimum and maximum values of the solution Method  Minimum  Maximum  Galerkin SUPG first order SUPG double asymptotic SUPG discontinuity capturing  -1.83421E-04 -3.78362E-02 -3.78362E-02 -8.77571E-04  1.33327 1.17226 1.17226 1.05377  The Galerkin method has the highest over-prediction in Table E1. Wiggles are most apparent with the method as show in figure E3a. The discontinuity capturing SUPG method (figure E3d) exhibits the least amount of wiggles although it is the most expensive in terms of computing time and effort.  177  (b) (a)  (d)  (c)  Figure E3. Three-dimensional plots of solutions from (a) Galerkin, (b) first order SUPG, (c) double asymptotic SUPG, and (d) discontinuity capturing SUPG  E.1.3 SUPG/Galerkin Comparison on Journal Bearings  Experimental data from Costa et al. [69] and Dowson [59] are simulated using the three methods. An estimate of the error is determined from (E.4)  178  where Tpred and Texp are the temperature predicted by the model and experimental temperatures respectively at each radial position on the bush. Table E2 and figure E4 show the results of the simulation.  Table E2 Comparison of SUPG and Galerkin predictions N , rpm  500 1000 T max T T max , T max , max ,Pred ,Exp. Error Pred. Exp. Error  Galerkin 41.0 SUPG 1 38.6 SUPG 2 38.9  38  0.078 51.6 0.016 45.4 0.024 45.5  43  0.199 0.057 0.058  T max 56.6 51.3 52.2  1500  2000  3000 4000 T max T max T max T max , T max T max , T max , Error Error Error ,Pred Error Pred. ,Exp. Pred. Exp. ,Exp. ,Exp. 0.098 59.0 51.5 0.004 52.9 0.014 52.9  0.110 67.9 53.1 0.003 59.2 0.003 59.5  0.152 59.0 0.005 0.010  76.6 65.3 65.8  68.0  0.126 0.040 0.032  Results indicate a steady increase in error for Galerkin method. The classic SUPG and double asymptotic methods give progressive better results with increasing speed. Although the latter method requires more computational effort in terms of time and speed, it appears there is no benefit in using it except for the highest speed of 4000 rev/min considered here. There is no doubt, however, of the benefit of better thermal predictions of SUPG especially for higher speed simulations. Considering that bearing speeds of 18000 rev/min are not uncommon, SUPG should be the method of choice in such cases. This, to the best of knowledge, represents the first such analysis in lubrication problems. More validation simulations of data in the literature follow in subsequent sections.  179  Performance of SUPG versus Galerkin 80 75 70  Temperature, C  65 60 55 50  Experimental  45  Classic SUPG  Galerkin SUPG: Doubly Assymptotic  40 35 30 0  500  1000  1500  2000  2500  3000  3500  4000  4500  Speed, rpm  Figure E4 SUPG and Galerkin predictions as a function of speed  E2  Slider bearing simulation  This example is from Arnell [66] (example 5.4.3, page 136) and analyzes three slider bearings using the infinitely long approximations (ILA) for maximum pressure pm, flow rate Qi, load capacity Wd, and frictional power loss PL. This example is from Arnell [66] (example 5.4.3, page 136) and analyzes three slider bearings using the infinitely long approximations (ILA) for maximum pressure pm, flow rate Qi, load capacity Wd, and frictional power loss PL. E2.1  Problem Statement  A thrust bearing of length L=0.2 m and width B=0.05 m operates with a minimum film thickness of 0.0001m and sliding velocity U=3 m/s, effective viscosity is µ=0.04 Pa.s,  180  and three tilt ratios K=(h2/h1-1) of 1.19, 2.0, and 0.50. The tilt ratio is defined for the three cases by keeping h1 constant at 0.0001 m and varying h2. Operating characteristics maximum pressure pm, flow rate Qi, load capacity Wd, and power loss PL are to be determined. Figure E5 shows the setup.  Stationary pad h2 y  h h1  z  U  x B  Runner  Figure E5 Slider geometry. E.2.2 Solution Procedure  Input data and results from ILA and SUPG FE are compared in table E3. The ILA expressions for maximum pressure pm, flow rate Qi, load capacity Wd, and power loss PL are [66] (E5) (E6)  ∆  181  1  (E7)  1  (E8) (E9)  Deviation values εr are defined as the difference between SUPG FE prediction and ILA divided by the SUPG figure.  Table E3 Prediction comparison, analytical vs SUPG FE Bearing Input U, m/s  μeff (Pa.s)  h1 (m)  Performance Parameters k 0.33  3  0.04 0.0001 0.46 0.67  p max (Pa) 150000 147879 153293 150914 120000 118201  QL, εr -0.014 -0.016 -0.015  m3 4.50E-05 4.5E-05 4.12E-05 4.12E-05 3.60E-05 3.65E-05  εr 0.010 0.001 0.012  Wd N 887.5 742.3 961.5 830.4 787.0 650.4  εr -0.196 -0.158 -0.210  PL, W 25.1 24.1 27.2 26.3 30.4 29.8  Δ T, εr -0.043 -0.031 -0.019  °C 0.34 0.32 0.4 0.39 0.34 -0.32  εr -0.043 -0.032 2.043  ILA SUPG ILA SUPG ILA SUPG  The values of εr on pressure, flow rate, and power loss are small (between 1 and 3%) which is very good. L/B ratio for this bearing is 4 which is not a long bearing, so ILA over-predicts the performance parameters.  E3  Journal Bearing Simulation  E3.1  Problem Statement  Arnell [66] (page 169): A hydrodynamic journal bearing of length 160 mm operates with a shaft 200 mm diameter which rotates at 1200 rev/min. The effective viscosity is assumed to be 40 cP (centipoises) and diametrical clearance is 0.18 mm. Determine the operating characteristics of the bearing at an eccentricity ratio of 0.7. E3.2  Solution Procedure  The problem is analyzed using SUPG FE and results compared to the analytical method used in Arnell. The method [68] uses tabulated numerical solutions for performance characteristics flow rate Qi, frictional power loss PL, and maximum pressure pm. Table E4 shows input parameters and predictions from the two methods.  182  E3.3 Simulation Results  The deviation between the two methods as defined in equation E4 above is small (εr<2.5%) for maximum pressure, flow rate and power loss. Temperature increase shows a deviation of 5.8%. The deviation between the two methods as defined in equation E4 above is small (εr<2.5%) for maximum pressure, flow rate and power loss. Temperature increase shows a deviation of 5.8%. Table E4 Input parameters and predictions. Bearing Input  E4  μ eff , Pa.s  0.04  k f , W/m.K C, m  0.13 9.00E-05  U, rev/s  20  L, m  0.16  D, m C p , J/kg.K  1880  ε  0.7  p max Pa  Performance Parameters QL PL  ΔT  2.64E+07 1.38E-04  8988  31.40 ILA  2.64E+07 1.40E-04  9193  0.0002  m 3 /s  0.0143  0.0223  29.68 SUPG -0.058  εr  0.2  Bearing data from Mitsui  Mitsui’s work [72] represents a landmark in journal bearing thermal predictions. The work presents bearing thermal performance as a function of speed Ns, clearance ratio C, type of lubricant, load capacity. Three different oils (transformer, no. 90 and no. 140) with different absolute viscosities are used. Experimental data and model predictions in figure E6 show the accuracy of the proposed model. Mitsui’s predictions on no. 90 and 140 oil are included for comparison. SUPG appears to give reasonable results in all cases. Of significance is the improved boundary conditions imposed in the groove/mixing zone. Conservation of mass flow (equation 2.24) is imposed, resulting in better predictions than  183  in the results of section 4.3 and figure 4.12. Model predictions are still least accurate in converging region of the bearing (0<θ<150°) similar to the pattern in figure 4.12. It appears Mitsui’s model under-predicts the same region.  Table E5. Input values for Mitsui bearing Input Parameter  Mitsui [72]  Bearing Radius R, m Bearing width L, m Bush outer radius R2(assumed), m Radial clearance at 20°C, C, m-6 Groove angle, ° Groove length, m-3 Groove Location, ° Specific heats cf, cb, cs Conductivities kf, kb, ks, kg, W/m K Convective heat transfer coefficient , bush/ambient H1, W/m2.K Convective heat transfer coefficient , bush/oil and shaft/oil interfaces H2, W/m2.K Inlet temperature Ti, K Ambient Temperature To, K Inlet lubricant density ρI, kg/m3 Inlet viscosities μi, Pa.s  0.05 0.07 0.1 78.5 10 60 0 1940, 380, 490 0.13, 250, 50, 0.025  Viscosity-temperature coeff. β, /°K Supply pressure ps (gauge), kPa  184  80  250 313 300 865 0.0736 (transf. oil) 0.0192 (No. 90 oil) 0.0469 (No. 140 oil) 0.0414 98  Predictions vs Exp. Temperatures 60  58 Mitsui:No.140 No.140  56  No.90 Transformer SUPG:Transf  54  SUPG:No.90  Temperature, C  SUPG: No.140 Mitsui:No.90  52  50  48  46  44  42  40 0  30  60  90  120  150  180  210  240  270  300  330  360  Angle, degrees  Figure E6. SUPG predictions versus experimental temperatures, Mitusi [72]. Included is Mitsui’s prediction for no.90 oil.  A summary of the maximum and average errors for SUPG and Mitsui’s model are shown in table E6. SUPG predictions are significantly better in each case.  Table E6 Maximum and average error for Mitsui’s data. Error Summary Transformer oil No. 90 oil SUPG Mitsui SUPG Mitsui Max Error Ave. error  0.017 0.008  0.039 0.017  0.014 0.026  185  0.030 0.064  No. 140 oil SUPG Mitsui 0.015 0.025  0.028 0.073  E5. SUPG performance in predicting Pierre’s data  The data and THD model presented by Pierre and Fillon [73] is analyzed and compared to SUPG predictions. The researchers present thermal and performance parameters for bearing input data presented in table E7.  Table E7 Input parameters for journal tested [73]  Input Parameter [unit] Diameter D, [m] Bearing length L, [m] Radial clearance C, [µm] Rotational speed Ns, [rev/min] Oil supply temperature Ti, [°C] Bearing load Wd, [N] Dynamic viscosity µì40°C at 40 °C (10-3 Pa.s) Dynamic viscosity µì100°C at 100 °C (10-3 Pa.s) Lubricant density ρ, [kg/m3] Lubricant specific heat Cp, [J/kg.K]  Value 0.1 0.08 123 3000 35 and 50 9000 30 4.05 860 2000  Experimental and predicted temperatures for the two speeds are shown in figure E7 (figure 2 in reference 73). Their THD model appears to be less accurate in the convergent and divergent zones.  186  Pierre's Model  70 T35Exp. T50Exp. Pierre@35C Pierre@50C  Temperature, C  65 60 55 50 45 40 35 0  60  120 180 240 Angle from inlet groove, deg  300  360  Figure E7. Experimental and THD temperature predictions as a function of angle from inlet groove [73].  The two cases in figure E7 are analyzed using SUPG and presented in figure E8. Included on figure E8(a) is Pierre’s prediction for reference. Table E7 gives a summary of the error magnitudes from the two methods. These show SUPG to give better results in both cases. Both models are off on the convergent region while SUPG gives more accurate predictions in the divergent zone (240<θ<360). Table E8 gives a summary of performance parameter predictions for Pierre’s data. L/D ratio is 1, input oil temperature Ti is 40°C (other input parameters as shown in table E6). The bearing characteristic number S. The error, as determined as before (equation E3) is small, showing the accuracy of SUPG.  187  SUPG vs Pierre Predictions, Ti=35 C 65  60 Exp. Data SUPG Pierre  Temperature, C  55  50  45  40  35 0  30  60  90  120 150 180 210 240 270 300 330 360 Angle from groove pos., degrees  (a) SUPG vs Pierre Predictions, Ti=50 C 70 65  Exp. data SUPG Pierre  Temperature, C  60 55 50 45 40 0  30  60  90 120 150 180 210 240 270 300 330 360 Angle from inlet groove, degrees  (b) Figure E8. Comparison of SUPG and Pierre’s THD model [73]. 188  Table E8. Error on thermal predictions, SUPG vs Pierre’s THD model  Max Error Ave. error  Error Summary T50 T35 SUPG Pierre SUPG Pierre 0.047 0.052 0.058 0.076 0.105 0.128 0.029 0.037  Table E9. Performance parameter predictions for Pierre’s data S=0.227 PL  Exp. Data 820  SUPG 819.5  Error 0.0006  pmax, MPa  2.3  2.26  0.0174  58 4.50E-05  59.5 4.60E-05  0.0250 0.0222  Tmax, C QL, m3/s  189  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0067254/manifest

Comment

Related Items