UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 ANALYSIS AND DESIGN OF PLANAR MULTIBODY SYSTEMS WITH REVOLUTE JOINT WEAR By SAAD M. MUKRAS A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2009 PAGE 2 2 2009 Saad M. Mukras PAGE 3 3 To my parents, Mohamed Mukras and Bauwa M ukras, to my siblings AbduRahman, Suleiman, and Mariam, to my wife Amina and finally to my son Talha PAGE 4 4 ACKNOWLEDGMENTS I expres s my humility and utmost gratitude to Allah for his blessings in my life. Verily no success would have been achieved without his gra ce and mercy. I would next like to thank my parents for their support in all as pects of my life. I owe them mu ch more than I can ever give back. Next I thank my wife for her unwavering support, her encouragement and for the patience that she has shown as I pursue my studies. I would like to acknowledge Dr NamHo Kim, my adviser, for the support that he has provided. Because of his advice and challenges, I have matured as a student and as a re searcher. I would like to thank my graduate committee members, Dr Tony L. Schmitz, Dr W.G. Sawyer, Dr B.J. Freg ly and Dr Jrg Peters for their assistance and guidance during my Ph.D. pursuit. I would also like to acknowledge the assistance that I have received from my colleagues, friends and member s of the university staff. Indeed, it would be negligent not mention the support that I have received from the members of the Masaajid in Gainesville who have enabled me to feel at home while away from home. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........8 LIST OF FIGURES.......................................................................................................................10 NOMENCLATURE......................................................................................................................13 ABSTRACT...................................................................................................................................14 CHAP TER 1 INTRODUCTION..................................................................................................................16 Motivation...............................................................................................................................16 Background.............................................................................................................................17 Objectives...............................................................................................................................20 2 DYNAMIC ANALYSIS OF RI GID MULTIBODY SYSTE MS.......................................... 22 Introduction................................................................................................................... ..........22 Kinematics Analysis............................................................................................................ ...23 Kinematic Constraints..................................................................................................... 24 Examples of Kinematic Constraints................................................................................ 26 Revolute joint...........................................................................................................26 Translational joint.................................................................................................... 27 Example: Kinematic Analysis of a SliderCrank Mechanism................................................ 28 Dynamic Analysis............................................................................................................... ....29 Direct Integration............................................................................................................. 31 Constraint Stabilization Method...................................................................................... 32 Generalized Coordinate Partitioning Method .................................................................. 33 Hybrid Constraint StabilizationGenera lized Coordinate Partitioning Method .............. 34 Modified Lagrangian Formulation.................................................................................. 34 Example: Dynamic Analysis of a SliderCrank Mechanism.................................................. 34 Summary and Discussion.......................................................................................................35 3 DYNAMICS OF MULTIBODY SYSTEMS WITH IMPERFECT REVOLUTE JOINTS ...................................................................................................................................44 Introduction................................................................................................................... ..........44 General Imperfect Revolute Joint........................................................................................... 46 ContactImpact Force Model...........................................................................................46 Modeling the General Imperfect Joint Model................................................................. 51 PAGE 6 6 Example: SliderCrank Mechanism with Joint Clearance Between the Crank and the Connecting Rod (General Im perfect Joint)............................................................ 54 Summary and Discussion.......................................................................................................56 4 WEARPREDICTION METHODOLOGY........................................................................... 62 Introduction................................................................................................................... ..........62 Wear Model..................................................................................................................... .......64 Wear Simulation Procedure....................................................................................................67 Computation of Contact Pressure....................................................................................67 Determining the W ear..................................................................................................... 68 Geometry Update Procedure........................................................................................... 69 Boundary Displacement..................................................................................................72 Reducing Computational Costs.............................................................................................. 73 Adaptive Extrapolation Scheme............................................................................................. 76 Experimental Validation........................................................................................................ .77 Summary and Discussion.......................................................................................................78 5 INTEGRATED MODEL: SYSTEM DYN AMICS AND W EAR PREDICTION................ 88 Introduction................................................................................................................... ..........88 Dynamic Analysis............................................................................................................... ....88 Wear Analysis.................................................................................................................. .......89 Demonstration of the Integration Process...............................................................................90 Summary and Conclusions.....................................................................................................92 6 INTERGRATED MODEL: SYSTEM DYNAMICS AND WEAR PR EDICTION USING THE ELASTIC FOUNDATION MODEL...............................................................97 Introduction................................................................................................................... ..........97 Elastic Foundation Model.......................................................................................................97 Analysis of Multibody Systems with Joint Wear Using the EFM........................................ 100 Modeling the NonIdeal Revol ute Joint Using the EFM ............................................... 100 Integrated Model: System Dynamics and W ear Prediction Using EFM....................... 103 Summary and Conclusion.....................................................................................................106 7 EXPERIMENTAL VALIDATION OF THE INTEGRATED MODELs ............................ 113 Introduction................................................................................................................... ........113 Experiments for Model Validation....................................................................................... 113 Summary and Conclusions...................................................................................................116 8 DESIGN OF A MULTIBODY SYSTEM FOR REDUCED JOINT W EAR MAINTENANCE COSTS.................................................................................................... 124 Introduction................................................................................................................... ........124 Design Example: Design of a Slider Crank for Reduced Maintenance Cost .......................125 Problem Definition........................................................................................................ 125 PAGE 7 7 Analysis of the SliderCrank w ith Mu ltiple Joints Wearing......................................... 126 Solution of Optimization Problem................................................................................. 127 Summary and Concluding Remarks.....................................................................................130 9 SUMMARY AND FUTURE WORK.................................................................................. 137 Summary and Discussion.....................................................................................................137 Future Work..........................................................................................................................140 LIST OF REFERENCES.............................................................................................................142 BIOGRAPHICAL SKETCH.......................................................................................................151 PAGE 8 8 LIST OF TABLES Table page 21 Dimension and mass parameter for slidercrank mechanism................................................. 37 31 Dimension and mass parameter for slidercrank mechanism................................................. 58 32 Material properties for the joint com ponents......................................................................... 58 33 Parameters for the force model............................................................................................ ...58 41 Wear test information for the pinpivot assembly .................................................................. 81 42 Simulation parameters for the pinpivot sim ulation test........................................................81 43 Comparison of results form th e sim ulation and Expt. wear tests........................................... 81 51 Dimension and mass properties of the slidecrank mechanism.............................................. 93 52 Properties of the pin and bushing...........................................................................................93 53 Test and simulation parameters............................................................................................ ..93 61 Dimension and mass parameter for slidercrank mechanism............................................... 107 62 Material properties for the joint com ponents....................................................................... 107 63 Wear simulation parameters................................................................................................ .107 64 Comparison of results from prediction based on FEM and EFM......................................... 107 71 Dimension and mass properties of the slidecrank mechanism............................................ 119 72 Properties of the pin and bushing.........................................................................................119 73 Test and simulation parameters............................................................................................ 119 74 Comparison of wear results for FE M and EFM m odels (21,400 crank cycles)................... 119 75 Comparison of wear results between test and FEM m odel (21,400 crank cycles)............... 119 76 Comparison of wear results between test and EFM m odel (21,400 crank cycles)............... 120 81 Dimension and mass parameter for slidercrank mechanism............................................... 132 82 Material properties for the joint com ponents....................................................................... 132 83 Parameter and design space specifications for the optimization.......................................... 132 PAGE 9 9 84 Solution of optimization problem (Eq. (81)).......................................................................132 85 Comparison of results between surrogate and highfidelity m odel...................................... 132 PAGE 10 10 LIST OF FIGURES Figure page 11 Examples of Multibody systems............................................................................................. 21 21 Revolute joint between bodies i and j .....................................................................................38 22 Translational joint between body i and j................................................................................ 38 23 A slidercrank mechanism.................................................................................................. ....39 24 Disassembled slidercrank mechanism................................................................................... 39 25 Results from the kinematic analysis (crank)...........................................................................40 26 Results from the kinematic analysis (slider)........................................................................... 41 27 Direct integration procedure fo r the Differential Algebraic Equation ................................... 42 28 Effect stabilization parameters ( and ) on the second kinematic constraint (2 K )...........42 29 Results from the dynamic analysis......................................................................................... 43 210 Joint reaction force between th e crank and the connecting rod............................................ 43 31 A revolute joint with clearance........................................................................................... ....59 32 KelvinVoigt viscoelastic model............................................................................................59 33 Contactforce model with hysteresi s dam ping (by Lankarani & Nikravesh)......................... 59 34 General imperfect joint................................................................................................... ........60 35 Penetration during contact between the pin and the bushing ................................................. 60 36 Slidercrank mechanis m s with joint clearance....................................................................... 61 37 Comparison of reaction between the ideal and the nonideal joints ....................................... 61 41 Pinpivot assembly........................................................................................................ .........82 42 Pinpivot finite element model............................................................................................ ...82 43 Wear simulation flow chart for the step update procedure .................................................. 83 44 A threenode contact element used to represent the contact surface ...................................... 83 45 Geometry updates process.................................................................................................. ....84 PAGE 11 11 46 Geometry updates process usin g the boundary displacement approach.................................84 47 Contact pressure distri bution on a pinpivot assem bly........................................................... 84 48 Extrapolation history for a pinpivot assem bly...................................................................... 85 49 Cumulative maximum wear on pin and pivot........................................................................ 85 410 Evolution of the contact pressure for nine interm ediate cycles............................................ 86 51 Integration of wear analysis into system dynamics analysis.................................................. 94 52 Slidercrank mechanism with a wearing jo int between the crank a nd the connecting rod .... 94 53 Initial joint reaction force for joint with clearance 0.0005mm............................................... 95 54 Locus of contact point C for a com plete crank cycle........................................................... 95 55 Wear on bushing as a function of the bus hing angular coordina te after 5,000 cycles ........... 96 61 Elastic foundation model.................................................................................................. ....108 62 Procedure to determine the contact pressure using the EFM............................................... 108 63 General imperfect joint................................................................................................... ......109 64 Penetration during contact between the pin and the bushing ............................................... 109 65 Slidercrank mechanis m s with joint clearance..................................................................... 109 66 Comparison of reaction force between the ideal joint m odel and the EFM......................... 110 67 Locus of contact point C for a com plete crank cycle......................................................... 110 68 Comparison of reaction force between ideal joint m odel EFM............................................ 111 69 Comparison of joint wear betw een the EFM and FEM after 5,000 cycles .......................... 112 71 Experimental slider crank mechanism.................................................................................. 121 72 Comparison of the initial jo int reaction force between th e two models and the Expt.......... 121 74 Comparison of the wear prediction between the m odels......................................................122 75 Comparison of the wear profile for the m odels and the experiment.................................... 123 81 Backhoe system with three main revolute joints.................................................................. 133 82 Slidercrank mechanism with two wearing joints wearing .................................................. 133 PAGE 12 12 83 Joint reaction force for Joint1 and Joint 2 ............................................................................ 134 84 Incremental sliding a ngle Joint 1 and Joint 2 ....................................................................... 134 85 Locus of the center of the contact region............................................................................. 135 86 DOE used to construct surrogate for cyc1 and cyc2..............................................................135 87 Response generated using the surrogate for cyc1..................................................................136 88 Response generated using the surrogate for cyc2..................................................................136 PAGE 13 13 NOMENCLATURE A Contact area E A Extrapolation factor Penetration re Coefficient of restitution E Youngs modulus WE Composite elastic modulus (Elastic foundation model) N F Normal force in the contact interface h Wear depth k Dimensioned wear coefficient K Nondimensional wear coefficient L Thickness of elastic layer (Elastic foundation model) Vector of Lagrange multipliers M Mass matrix p Contact pressure q Position vector AQ Vector of applied loads s Sliding distance t Time Poisons ratio Constraint vector q Jacobian matrix Crank velocity PAGE 14 14 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ANALYSIS AND DESIGN OF PLANAR MULTIBODY SYSTEMS WITH REVOLUTE JOINT WEAR By Saad M. Mukras August 2009 Chair: Nam Ho Kim Major: Mechanical Engineering Wear prediction on the compone nts of a mechanical system without considering the system as a whole will, in most cases, lead to in accurate predictions. This is because the wear is directly affected by the system dynamics whic h evolves simultaneously with the wear. In addition, the contact condition (regions of contact for the wearing bodies) also depends on the system dynamics and, in most cases, can only be determined in a multibody dynamics framework. In this work, a procedure to analyze planar multibody systems in which wear is present at one or more revolute joints is presented. The analysis involves modeling multibody systems with revolute joints that consist of cl earance. Wear can then be incorporated into the system dynamic analysis by allowing the size and shape of the clearance to evolve as dictated by wear. An iterative wear prediction proce dure based on the Archards wear model is used to compute the wear as a function of the evolving dynamics a nd tribological data. In this framework, two procedures for the analysis of planar multibody sy stems with joint wear are developed. In the first procedure contact force at the concerned jo int is determined using a contact force law and the wear prediction is based on the finite elemen t method. In the second procedure, contact force determination and the wear prediction are based on th e elastic foundation model. PAGE 15 15 The two procedures are validated by comparing the wear predictions with wear on an experimental slidercrank mechanism. The experime ntal slidercrank is also used as a reference to assess the performance of the two models. It turns out that the procedure based on the finite element method provides reasonably accurate pr edictions for both wear profile and wear volume/mass whereas the procedure based on th e elastic foundation model provides reasonably accurate estimates on the wear volume/mass, is co mputationally faster but provides progressively poor estimates on the wear profile. Finally an example is presented to illustrate an application of the procedures. In the example, a slidecrank mechanism is designed so that the maximum allowable wear depth at two joints occurs simultaneously. PAGE 16 16 CHAPTER 1 INTRODUCTION Motivation In this work, a study of the analysis of plan ar multibody systems with wear at a single or multiple joints is presented. This analysis will involve the dynamic analysis of multibody system coupled with a wear prediction procedure to estimate the wear that the wearing joints. The motivation behind this work is based on the need to understand the behavior of such systems under the influence of wear. Informed and theref ore improved design of these systems is thus enabled. Figure 11 shows two examples of multibody /m echanical systems namely a backhoe system and a flat engine. The two systems have different functionalitie s but are, otherwise, similar in the sense that they are both made of multiple components connected by joints. Due to the contact and relative motion of the joint components it is inevitable that wear will occur at the joints of both systems. Despite this fact, knowledge of how wear affects the system dynamics and how the system dynamics affect the wear ca n empower the designer to develop better design. For instance, in the flat engine shown in Figure 11, it is known from the analysis of a sim ilar system (in this work) that the revolute join t between the crank and the connecting rod will experience the most wear. This is especially true if the reaction force at this joint and that of other joints are similar (which is the case sin ce the resultant reaction fo rce on the piston due to compression in the cylinder will dom inate inertial forces). One po ssible way to reduce the wear is to increase the joint diameter; this will ha ve the effect of reducing the maximum contact pressure and thus reducing the maximum wear depth. On the cont rary increasing the diameter will have the opposing effect or increasing the wear since the incremental sliding distance is proportional to the diameter. Thus the decision to increase or decrease the diameter will depend PAGE 17 17 on which of the two quantities (the contact pressure or the slidi ng distance) is more dominant. Arriving at this decision will require an analysis that couples both the system dynamics and wear prediction which will be the discussion of this work. Another example can be drawn from the backhoe system shown in Figure 11. Being able to perform multibody analysis on this system while accounting for wear at the joints, can allow the engineers to correctly predict wear at the joints and thereafter issue appropriate warranties or prepare appropriate maintenance schedules. Prediction of wear at the joints is however not trivial. Wear prediction at component level ma y be inaccurate since these procedures do not account for how the overall system will affect a particular joint. For instance, it will not be a simple task to determine the regions in a joint that will be in contact and thus wear without considering the entire system. In addition wear th at occurs at the joint may further affect the subsequent regions of contact. On ce again it is clear that to determine this information it is necessary to perform an analysis in which the system dynamics is coupled with a wear prediction procedure. Background As was m entioned earlier, the analysis of mu ltibody systems with joint wear will be the focus of this work. While this work is not uni que, more research is st ill required to properly address this area. Most of the work related to th is topic has been aimed at modeling the effect of joint clearance on multibody systems. Early studies, in this area, focused on simple models to obtain insight into the behavior of systems with joint clearances. Dubousky et al [1,2] developed a contact impact pair model to study the elastic joint with clearance. In their model, joint elasticity and damping were modeled via springs a nd a viscous coefficient. The study served to represent the complexity of the joint clearance using a simple model. Dubowsky and Gardner [3] later extended this work to include flexible m echanism as well as multiple clearance connections. PAGE 18 18 Earles et al [4] proposed to model joint cl earance using a massless rigid link whose length was equal to the clearance size. The components of the joint were thus assumed to be in contact at all times. Wu et al [5] later used the model to predict contact loss between the joint components for planar mechanism. The concept of a massless link was also used by Furuhashi et al [69] to study the dynamics of a fourb ar linkage with clearance. Once again the joint components were assumed to be in contact at all times. In later year, more complex models were de veloped to study the e ffect of clearance on system dynamics. Farahanchi et al [10] m odeled joint clearance by considering three configurations of the joint components i.e. 1) fr eeflight motion; when the components are not in contact 2) impact condition; when the component s establish contact and 3) sliding condition; when the components are in cont act and relative motion. In their model the reaction force at the joint (when the joint component s are in contact) was determ ined assuming no clearance was present. They used a slidercrank mechanism to demonstrate the procedure and studied the effect of clearance size, friction, crank speed and impact parameters. Rhee et al [11] also used the three modes of motion to model the joint clearance. Th ey used an approach similar to that of Farahanchi [10] to determine the reaction fo rce during the sliding motion. They studied the response of a fourbar mechanism w ith a revolute jo int clearance. The three mode approach was also used by Khulie f et al [12] to model the clearance at the joint. In their approach, termed as discon tinuous method, the analysis (integration of the equations of motion) is divided into two pa rts namely; pre and postcollision during which momentum balance is performed to determine th e postcollision velocities. Velocities in the analysis are then updated a nd the analysis is resumed. PAGE 19 19 Ravn [13], once again, utilized the threemode a pproach to model the joint clear. However, in his approach, the reaction for ce during the impact and the slid ing mode is computed using a contact force model. The analysis in this case has been termed continuous since integration of the equations of motion is not halted as in the case of the discontinuous method. A number of researchers [1418] have since used this techni que to model as well as study the effect of clearance in the joints of multibody systems. The literature that has been presented in the preceding paragraphs was primarily directed at modeling multibody systems with joint clearance th at remains constant. The concern of this work, however, is in studying the case in which th e clearance changes in shape and size due to wear. The majority of wear pred iction procedures that are avai lable generally take into account the changes in the contact pressu re distribution due to the geometry evolution caused by wear. This is however at a component level rather than at the system level wh ere the effect of the system dynamics is accounted for in the wear prediction. Recently, Flores [19] presented a procedure to model multibody system s with joint wear. The procedur e involves the analysis of multibody system in the manner sim ilar to the threemode appro ach. Based on the contact force at the concerned joint, wear is estimated for the joint. While this work is in line with the basic theme of this dissertation details of the work differ substantiall y. These include: 1) the manner in which the system dynamics is integrated with th e wear prediction, 2) co mputation of quantities required for wear analysis such as i) calculati on of the sliding distance, ii) prediction of the contact surface, and iii) calculation of the cont act pressure, 3) the procedure in which the geometry is updated to reflect the wear, 4) the procedure in which the analysis is accelerated in order to reduce computational expense, 5) the manner in which the contact force is computed, PAGE 20 20 and 6) experimental validation of the pro cedures. These differen ces have significant consequences in the performance of the procedures presented. Objectives In line with the m otivation, f our objectives to this work are listed as follows: To develop a procedure that can be used to analyze multibody systems with joint wear. To demonstrate the need of a multibody framew ork in predicting wear at the joints of mechanical systems. To validate the procedure through experiments. To present an example that will illustrate the use of the procedure as a design tool. Since the subject of multibody system analysis is extensive, this work will be limited to analysis of planar systems. Furthermore the components in the multibody systems will be assumed to be rigid. While in real system wear is expected to occur in all components that are in contact and in relative motion, this work will only consider wear at the revolute joints. It is, however, possible to extend this work to cover wear in other types of joints. PAGE 21 21 Figure 11. Examples of Multibody systems. Revolute joints Flat Engine Piston Connectingrod Crank Revolute joints Backhoe System PAGE 22 22 CHAPTER 2 DYNAMIC ANALYSIS OF RI GID MULTIBODY S YSTEMS This Chapter presents a discussion on the dynamic analysis of multibody systems. Formulation and solution techniques of the e quations of motion for in plane, rigid, multibody systems are presented. The intent of the information in this Chapter is to acquaint the reader with the subject of multibody dynamics and lay the fou ndation for the subsequent topic of analyzing multibody systems with joint wear. A more in de pth discussion of the su bject can be accessed from the literature [2024]. Introduction A m ultibody system refers to a system consisting of a set of interconnected rigid or flexible bodies that undergo large displacements and rotations. Many mechanical systems such as vehicles, robots, etc are in fact multibody systems. They consist of several interconnected bodies that may be loaded at various locations and w hose motions are restricted to achieve a desired function. The analysis of such systems to determin e their response to applied loads is termed as dynamic analysis of multibody systems. Depending on the degree of complexity of the system being modeled, the formulations of the equations of motion can be quite simple or complicated. The most general example of such a system woul d consist of a set of rigid and/or flexible components subjected to various lo ading conditions (forces contact force, torque etc) and whose motions are spatial. In reality, however, appropriate assumptions can be made so as to manage the complexity of the analysis for a corres ponding system. For inst ance, a systems whose motions are primarily confined to a single plane may be modeled as plan ar instead spatial and thus substantially reducing the modeling efforts. Another example would be to model a system with rigid bodies instead of flexible bodies if the deformations of the bodies are reasonable small. PAGE 23 23 The two simplifications previously men tioned are probably the most important consideration when analyzing a multibody system. In this work only systems with components whose motions are restricted within a plane are c onsidered. In this Chapter, dynamic analysis of systems with rigid bodies (referred to as rigid multibody dynamics) will be discussed. As was mentioned earlier, dynamic analysis of a multibody system is usually conducted in order to determine the response of a system to a pplied loads. Thus the joint reactions forces or the motion (position, velocity and acceleration) of the components of the system are typical quantities of interest. In cases in which only the motion of the system is require, it may be possible to simply perform a kinematic analysis. A discussion of dynamics analysis will thus be preceded by a brief outline of the kinematic analysis of a mechanical system. Kinematics Analysis Kinem atic analysis involves estimating the position, velocity a nd acceleration of the components of a multibody system without regard to the forces that produce the motion. The analysis efforts are therefore less than that of the dynamic analysis. The position and orientation of all bodies in a multibody system can uniquely be specified by a set of generalized coordinate s which may vary in time if the system is in motion. In this text, the generalized coordinates are expressed as a column vector, 12,,,T ncqqq q where nc is the number of generalized coordinates. Specif ication of the system components for a planar system is achieved by fixing an xy reference frame on each body in the system. Body i can then be located by specifying the global coordina tes of the origin of th e corresponding reference frame ,T i i x y r and the relative angle between the re ference frame and the global frame i The location of a body in the system is thus specified as ,,T i ixy q The generalized PAGE 24 24 coordinates can then be expressed as 12 12,,,,,,,,,,,T T TTT TTT nb nbxyxyxy qqqq where nb is the number of bodies in the multibody system It is thus clear that a planar system with nb bodies will have nc = 3nb generalized coordinates. Kinematic Constraints The bodies in a m ultibody system are interconne cted by joints which impose condition on the relative motion of the bodies. Consequently, the generalized coordinates are usually not independent. When these conditi ons are expressed as algebraic equations in terms of the generalized coordinates, they are referred to as holonomic kinematic constraints. The kinematic constraints for a system with nh constraints can be expressed as 12,,,T KKKK nh qqqq0. (21) Furthermore if the constraints depend on time t they are referred to as timedependent constraints and are expressed as ,Kt q0 (22) If t does not appear explicitly in Eq. (22) then the constraints are ca lled stationary constraints. Other general constraints that c ontain inequalities and/or depe nd on system coordinates as well as velocities are said to be nonholonomic constr aints. In this work only holonomic kinematic constraints are considered and are re ferred to simply as constraints. If the number of generalized coordinates exceed the number of constraints (ncnh ), then the constraint equations cannot be solv ed to uniquely determine the position ( q ) of the components. The system in this case is said to have ncnh degreesoffreedom ( DOF ). In order to determine the motion of the system one may specify a total of DOFncnh additional PAGE 25 25 constraints. These additional cons traints are termed as driving c onstraints and are shown in Eq. (23). ,Dt q0 (23) If the kinematic as well as the driving constraint s are independent the complete set of constraints (shown in Eq. (24)) will consist of nc(wherencnhDOF ) independent constraint equations which can simultaneously be solved to uniquely determine the position (tq) of the components of the system at any time. ,K Dt t t q q0 q (24) Alternatively, instead of specifying the driving constraints, one may specify appropriate loads on the system. In this case, tq is the solution of a set of differential equations (the differential equations of motion). This will be th e subject of the dynamic analysis. The expression for the velocity q as well as the acceleration q of the system components can be obtained by differentiating the constraint vector (Eq. (24)). Using the chain rule expression (24) is differentiated once with respect to tim e to yield the following velocity equation, t q q (25) where qis the Jacobian of the constraint vector. The velocity q can then be determined if the Jacobian is nons ingular. Equation (25) can further be differentia ted to yield the expression for the acceleration. The res ult is shown in Eq. (26). Similar to the velocity equation, if the Jacobian is nonsingu lar the acceleration can also be determined. PAGE 26 26 Determination of the position, velocity and a cceleration completes the kinematic analysis. It is worth noting that instead of supplying the dr iving constraints, approp riate driving forces or torques can be applied. In this case a kinematic analysis cannot be performed (since the numbers of unknowns exceed the number constraint equati on), and instead a dynamic analysis will have to be performed. 2 2ttt ttt qqq q qqq q qq q q 0 q qq q (26) Examples of Kinematic Constraints Various types of kinematic constraints can be formulated in order to achieve desired relative motion between components in a multi body system. In what follows, two kinds of constraints, the revolute joint and translational join t, will be presented in order to demonstrate the formulation procedure. The two constraints have been specifically selected as they will also be used in the later sections of this work. Revolute joint. The revolute joint imposes a constraint on the relative translation between two bodies i and j However, it permits relative rota tion between the bodies about a point P that is common to both bodies. An illustration of the revolute joint is shown in Figure 21. The constraint equation that descri b es the revolute joint can be fo rmulated by ensuring that point P on body i and point P on body j coincide at all times. This constr aint can be expressed in vector form as (,)rij iiijjj rAsrAs0, (27) where ir and jr are position vectors in the global coordinate system ( xy ) that describe the location of the origin of the bodyfixed coordinates ( xiyi and xjyj). The vectors isand js are PAGE 27 27 position vectors in the bodyfixed coordinate systems that locate the point P iA and jA are transformation matrices that transform vectors in bodyfixed coordinates sy stems into vectors in the global system. Translational joint. The translational joint, contrary to the revolute joint, allows relative translation between two bodies but prevents relative rotation be tween the bodies. For instance, Figure 22 shows two bodies, i and j which are connected by a translation joint. Body i trans lates relative to body j along the axis AB Relative rotation between th e two bodies is however constrained. With reference to Figure 22, this joint can be mode led by requiring that the vector ih on body i remain collinear to vector jh on body j The two will be collinear if ih is parallel to jh and ih is parallel to ijl, where ijl is a vector that c onnects a point on vector ih and another on vector ih. The constraint can be expressed mathematically as (,) iij tij ij hl 0 hh (28) It should be noted ijl in Eq. (28) can assume a value of zero when the two points in vector ih and ih coincide. This occurrence essentia lly means that the two vectors, ih and jh, touch each other and the second constraint (ij hh0) will ensure that the vectors are parallel. The two joints discussed above ar e simple but probably the most widely used constraints in planar multibody systems. The inte rested reader is referred the works of Nikravesh [20] and Haug [21], for a comprehensive disc ussion on constraint formulation. PAGE 28 28 Example: Kinematic Analysis of a SliderCrank Mechanism In order to demonstrate the procedure of cons traint formulation and kinematic analysis, a study will be done on a slidercrank mechanism. The slidercrank is se lected for this task since it is a simple mechanism with all the relevant features necessary for the demonstration. Furthermore the slidercrank will be used to facilitate the central idea in this work. The slidercrank mechanism, shown in Figure 23, consists of three rigid bodies (crank, connecting rod, and slider). In Figure 24, the disassembled components of the mechanism are shown in the global ax is. Each component can tran slate and rotate in the plane. The mechanism is modeled by imposing constraints on the motion of the components as described in the previous section. The constraints corresponding to the slidercrank mechanism, shown in Figure 24, consist of nine nonlinear simultaneous equations expressed explicitly as 11 111 21122 21122 31122 31122 3 3 1cos 0 sin 0 2coscos 0 2sinsin 0 2cos2cos 0 2sin2sin 0 0 0 0 xl yl xll yll xll yll y t (29) The first two constraints in Eq. (29) confine point P1 (on the crank) to the origin and describe a revolute joint between the origin/groun d and the crank. The next two constraints also describe a revolute joint between the crank and the connecting rod. They ensure that points P2 (on the crank) and P3 (on the connecting rod) coincide at all times. This condition is synonymous to an ideal joint and later will be relaxed when modeling the imperfect joint. The fifth and sixth constraints in Eq. (29) represent a perfect revolute join t be tween the conn ecting rod and the PAGE 29 29 slider. The next two constraints ensure that the slider remain on the xaxis w ithout rotation. These two constraints represent the translational joint. The final constraint is the driving constraint which specifies the motion of the crank. For the current case a constant angular velocity is imposed on the crank. It can be seen from the set of simultaneous equations above, that the number of equations is exactly equal to th e number of unknowns. The unknowns are the locations of the center of masses of the slidercrank components, denoted as 111222333,,,,,,,,Txyxyxy q. Kinematic analysis is performed by simultaneously solving Eq. (29) to determine q and solving Eq. (25) and (26) to determine the velocities q and accelerations q of the system. A kinematic analysis is performed for the slidercrank mechanism with crank velocity of rad/sec and dimension parameters as shown in Table 21. Figure 25 shows the position velocity and acceleration at the crank m a ss of center from the kinemati c analysis. The corresponding results for the slider are shown in Figure 26. Since the motion of th e slider is rest ricted along the global x axis (i.e. translation along the xaxis), the vertical position, velocity as well as acceleration of its mass center s hould all be zero. This is confirmed from the plots in Figure 26. Dynamic Analysis As was mentioned earlier, dynamic analysis is the study of motion and the forces that cause it. With reference to this definition, our primary interest is to develop and solve the equation of motion for the concerned system. Tw o approaches can be used to formulate the equations of motion namely: the NewtonEuler form ulation and the Lagrangian formulation. The NewtonEuler formulation involves the application of the law of motion which involves applied loads and reactions, whereas in the Lagrangian fo rmulation the system dynamics is described in PAGE 30 30 terms of work and energy. Both approaches lead to equivalent results. A more comprehensive discussion of the approaches can be found in th e works of Shabana [22] and Schiehlen [24]. The equations of motion for a constrained ri gid multibody system are stated here without derivation due to the length and complexity of thei r derivations. The interested reader is referred to the work of Haug [21] for a detailed deri vation of these equations using the Lagrangian approach. The equation can be expressed as TAqMq+ Q (210) where M is the mass matrix consisting of masses and moments of inertia for the system components. q and T q are the acceleration vectors and J acobian of the constraints (as previously defined) respectively. is a vector of Lagrange multipliers and AQ is a vector of applied load. It is worth noting th at in this work the bodyfixed coordinate system is selected to be at the center of mass (COM) of the corresponding body which significantly simplifies the general form of the equations of motion. In pa rticular, the mass matrix becomes diagonal. The product of the Jacobian and the ve ctor of Lagrange multipliers (Tq) is in fact the vector of reaction forces. This is the second term on the left hand side of Eq. (210). Thus for an unconstra ined system the vector or Lagrange multip liers are zero and this term would disappear. Equations (210) and (26) can be combined to result in a mixed system of differential algebraic equations also known as the differential algebraic equation of motion (DAE). The equation s are expressed as T A q qM q Q 0 (211) PAGE 31 31 where 2tttqqq q q qq q For a meaningful system th e coefficient matrix of Eq. (211) must be nonsingular. This is guaranteed if the m ass matrix M is positive definite and the Jacobian matrix q is full row rank (or constraints are independent). The process of determining the system dynamics involves solving the DAE. In addition to the DAE Eq. (211) both the position and velocity equation s must be satisfied, that is t q0, (212) t t q q q 0 q (213) In the next section a procedure termed as direct integration will be used to demonstrate how to solve the DAE (Eq. (211)). Direct Integration For a proper set of initial conditions (initial po sition and velocity conditions), the system of equations (Eq.(211)) can simultaneously be solved for a unique solution. The system of equations for a particular time instant consist of a set of linear equations that can easily be solved to determine the acceleration vector q and vector of Lagrange multipliers (which give the joint reactions). For that time instant, the velocity q and position q vectors can be determined by integrating the acceleration vector q using methods such as R ungeKutta, AdamsBashforth, and AdamsMoulton predictorcorrector methods [25]. The process is then repeated for the next time instant with the new values of position veloc ity and acceleration. This procedure is referred to as direct integration and is illustrated in Figure 27. Although the direct integration procedure is straight forward, th ere is no guarantee that, after several iterations, the positi on a nd velocity equations will be satisfied. This is especially PAGE 32 32 true because not all the components of q and q are independent. In order to combat this problem a number of algorithms have been deve loped. Some of these include the Constraint Stabilization Method, Generalized Coordinate Partitioning Method, Hybrid Constraint StabilizationGeneralized Coordinate Partitio ning Method, and the Modified Lagrangian Formulation. A brief discussion of these techniques will follow in the next sections. Constraint Stabilization Method The Constraint Stabilization Method, due to Baumgarte [26], is probably the most attractive procedure due to its simplicity both in concept and implementation. Baumgarte noted that the differential equations 0 can be unstable. Baumgarte, however, observed that if the acceleration equation, repeated here as q q 0 (214) was modified to 22 0 (215) such that 0 and 0 it became stable. It is clear from Eq. (215) that since the acceleration equation is zero ( 0) both the position and velocity equations are approximately zero i.e. 0. Equation (215) can explicitly be written as 22tqq q q (216) Both Eq. (214) and (215) are equivalent analytically howev er with r egards to stability the term 22 in Eq. (215) serves as an error control so that both the position and velocity constra ints will eventually be satisfied. PAGE 33 33 The DAE can then be modified by replacing th e acceleration equations with the modified equivalent. The resulting DAE is shown in Eq. (217). The system dynamics can then be determ ined by a similar procedure as outlin ed for the direct integration method. 22A T t q q qQ M q 0 q (217) The choices for and are important in the solution process, however, no generally accepted procedure appears to have been adopted in the selection of these parameters. Baumgarte [26] used a trial and error procedure through numerical experiments to determine reasonable values for the parameters. Lin et al [27,28] suggested that the appr opriate values of the parameters were dependent on the type of integration method empl oyed. They discussed ways of determining appropriate values of and for the RungeKutta method [27], the AdamsBashforth and the AdamsMoulton predictorcorre ctor methods [28]. Nikravesh [20] suggested that for most practical problems, values between 1 and 10 will suffice for both parameters. It was also suggested by Wittenbur g [23] that selecting will stabilize the response more quickly. Generalized Coordinate Partitioning Method The Generalized Coordinate Partitioning Method is a more rigorous technique that has been found to have good error control. The cons traints (position and velocity equations) are satisfied to the precision set by the user. The method involves partitioning the generalized coordinates into dependent and independent variables. Equation (211) can then be solved to determ ine the independent accelerations (together with the dependent accelerations) and the Lagrange multipliers The independent accelerations are in tegrated for independent velocities and positions. The kinematic constraint can then be solved to determine the dependent positions (using techniques such as the NewtonRaphson method). Finally the dependent velocity is PAGE 34 34 determined by solving the partiti oned velocity equation. More deta ils of this technique can be found in the works of Haug [21] and Wehage et al [29]. Hybrid Constraint StabilizationGeneralized Coordinate Partitioning Method As the name suggests, the Hybrid Constraint StabilizationGeneralized Coordinate Partitioning Method is a hybrid between the two techniques. Th e idea behind the method is to derive benefits from the two techniques na mely: to maintain good error control of the partitioning method and to retain the computational speeds of the constraint stabilization method. A complete study of the method is ava ilable in the works of Park [30]. Modified Lagrangian Formulation The approach used in the modified Lagr angian formulation method to remedy the constraint violations is quite diffe rent from the three procedures ju st described. In this procedure the constraint violations are ca tered for using a penalty appro ach. The constraint condition is inserted into the Lagrange equations by means of a penalty formulation rather than appending them as was the case in the nonmodified formulation (Eq. (211)). The benefit of this technique is that no new unknowns are introduced which is c ontrary to the regular Lagrange formulation where the new unknowns are the Lagrange multipliers. The details of the formulations have been presented by Bayo et al [31]. Four procedures to solve th e equations of motions have been presented. All these techniques have been shown to produce acceptable re sults. In this work, however, the Constraint Stabilization Method has been adopted because of its simplicity in implementation. Example: Dynamic Analysis of a SliderCrank Mechanism A slidercrank mechanism, which was introduced earlier, will be used to demonstrate the dynamic analysis procedure. A diagram of the slidercrank mechanism is shown in Figure 23. PAGE 35 35 The kinematic constraint equation for the slid ercrank are identical to those used in the previous example (Eq. (29)). In this case, the crank is cons trained (driving constraint) to rotate at a constant angular velocity of rad/sec. The dim ensions and mass properties for the mechanism are shown in Table 21. The stabilization parameters u sed for this analysis were 250 The differential equations of motion can then be assembled as described in Eq. (217). For this case, the only applied load is the g ravitational force. The equations of motion may be solved through the procedure outlined in Figure 27 in conjunction with the c onstraint stabilization method. In Figure 28, the value of the second kinematic constraint (constrai nt corresponding to the joint b etween the crank and the ground) is plotted for differe nt constraint stabilization parameters. It is clear that wh en no constraint stabilization para meters are used (i.e. direct integration method) the constraint is violated and the violatio n increases as the integration progresses. A value of 250 was found to provide reasonably stability to the constraint. The motion (position, velocity and accelerat ion) of the center of mass of the crank, obtained from the dynamic analysis, is plotted in Figure 29. The results shown in the figure show close resem bles to those obtained from the kinematic analysis in Figure 25. In addition to the m otion, the joint reaction forces at the various revolute joints can be determined. As was mentioned earlier, this informa tion can not be obtained by performing a kinematic analysis. The reaction force at the revolute joint between the crank and connecting rod is shown in Figure 210. Summary and Discussion In this Chapter the procedure to perform kinematic and dynamic analysis for a planar rigid multibody system was presented. The kinematic anal ysis involves simultaneously solving a set PAGE 36 36 of nonlinear kinematic constraint s equations. The analysis provid es information regarding the motion of the bodies in the system ; i.e., the positions, ve locities and accelerations of the bodies. In order to determine the reaction forces, a dynami c analysis is required. In the dynamic analysis, a differential algebraic equation of motion needs to be assembled. The equations of motion can once again be solved to determine the motion as well as reaction forces of the system. In the solution process it is necessary to integrate the accelerations in order to determine both the velocity and the position. It has, however, been no ted that a direct integration procedure will lead to violations in the kinematic constraints. Various techniques, such as the c onstraint stabilization, generalized coordinate partitioning and the modi fied lagrangian formulations, have been proposed to mitigate or eliminate the constraint vi olation problem. In this work the constraint stabilization method was used. The kinematic and dynamic analysis procedures were demonstrated using the s lidercrank mechanism. PAGE 37 37 Table 21. Dimension and mass pa rameter for slidercrank mechanism Length (m) Mass (kg) Moment of inertia (kgm2) crank 1.00 10.00 45.00 connecting rod 1.75 15.00 35.00 slider 30.00  PAGE 38 38 Figure 21. Revolute joint between bodies i and j. The joint imposes a constraint on the relative translation between the bodies. Figure 22. Translational join t between body i and j. The joint imposes a constraint on the relative rotation between the bodies. y x P yi x i i y j x j j r i r j si s j y x hi P i ih P j hi li j j i yi xi yi xi B A PAGE 39 39 Figure 23. A slidercrank mechanism Figure 24. Disassembled slidercrank mechanism Connectin g rod Cran k Slider Ground Ground y3 x3 x1 y1 y2 x2 y x y x P1 y1 x1 P2 P3 y2 x2 P4 y3 x3 P5 Crank Connecting rod Slider mass length crank m1 l1 connecting rod m2 l2 slider m3 l3 PAGE 40 40 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 1.5 1 0.5 0 0.5 1 1.5 2 Time (sec)Position (m)Position of Mass Center for the Crank x1 y1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 4 3 2 1 0 1 2 3 4 5 6 Time (sec)Velocity (m/s)Velocity of Mass Center for the Crank dx1/dt dy1/dt 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 15 10 5 0 5 10 15 20 Time (sec)Acceleration (m/s2)Acceleration of Mass Center for the Crank d2x1/dt2 d2y1/dt2 Figure 25. Results from the ki nematic analysis. Position velocity and acceleration of the crank mass center. PAGE 41 41 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 1 0 1 2 3 4 5 6 7 8 Time (sec)Position (m)Position of Mass Center for the Slider x1 y1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 8 6 4 2 0 2 4 6 8 10 12 14 Time (sec)Velocity (m/s)Velocity of Mass Center for the Slider dx1/dt dy1/dt 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 40 30 20 10 0 10 20 30 Time (sec)Acceleration (m/s2)Acceleration of Mass Center for the Slider d2x1/dt2 d2y1/dt2 Figure 26. Results from the ki nematic analysis Position velocity and acceleration of the slider mass center. PAGE 42 42 Figure 27. Direct integr ation procedure for the Differential Algebraic Equation 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.06 0.05 0.04 0.03 0.02 0.01 0 0.01 Time (sec)2nd ConstraintConstraint =0 =0 =5 =5 =100 =100 Figure 28. Effect stab ilization parameters ( and ) on the second kinematic constraint (2 K ). Initial conditions 0t,0q, 0q Assemble M, q, AQ, Solve linear DAE T A q qM q Q 0 for q Assemble vel. & acc. Vector t q g q Integrate for pos. and vel. integratett t q gg q Time increment ttt endtt PAGE 43 43 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 1.5 1 0.5 0 0.5 1 1.5 2 time (sec)Position (m)Position of Mass Center for the Crank x1 y1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 4 3 2 1 0 1 2 3 4 5 6 Time (sec)Velocity (m/s)Velocity of Mass Center for the Crank dx1/dt dy1/dt 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 15 10 5 0 5 10 15 20 time (sec)Acceleration (m/s2)Acceleration of Mass Center for the Crank d2x1/dt2 d2y1/dt2 Figure 29. Results from the dynamic analysis. Position velocity and acceleration of the crank mass center. 0 1 2 3 4 5 6 1500 1000 500 0 500 1000 1500 2000 2500 Crank Position (rad)Joint Reaction Force (N)Joint Reaction Force (crank and connecting rod) Fx Fy Figure 210. Joint reaction force between the crank and the connecting rod. PAGE 44 44 CHAPTER 3 DYNAMICS OF MULTIBODY SYSTEMS WITH IMPERFECT REVOLUTE JOINTS In this Chap ter a discussion on how to m odel the dynamics of multibody systems with imperfect revolute joints will be presented. In pa rticular constraints will be developed to model the imperfect joint. Introduction In the previous Chapter a discussion on th e procedure to model multibody systems was presented, based on kinematic constraint formula tions for revolute and translational joints. The formulations assume that clearances at these jo ints are nonexistent. However in real mechanical systems, joint clearances do in fact exist. The cl earances are due to manu facturing constraints as well as the occurrence of wear. In the latter case the two com ponents of the joint, are in contact and in relative motion and as a result wear occurs so that the clearance size increases with time. Joint clearances have been noted to affect th e performance and servic e life of mechanical systems. This may be attributed to the increas ed vibration, excessive wear and dynamic force amplification as discussed by D ubowsky [32]. The clearances may also cause uncertainty in the orientations and positions of mechanisms that are designed to have pr ecise motions such as manipulators [33]. Due to the significance of the problem, numerous studies have been conducted with the goal of understa nding the response of these syst ems in the presence of joint clearances [111,1318,3241]. These studies have evolved from the analysis of less complex planar multibody systems [111,1318,32,33,35,36,38,40,41] to more complex spatial systems [34,37,39] as well as from rigid multibody analysis [49,14,1618,32,33,38,40] to flexible multibody analysis [3,34,36,37,42]. The studies have demonstrated that the presence of clearances alter the response of the system appreciably. PAGE 45 45 Figure 31 shows a diagram of a revolute joint w ith a join t clearance. The joint consists of a pin and a bushing (also referred to as a journa l and a bearing, respectively). The pin and the bushing are attached to either of the two bodies that share the re volute joint. This joint differs from the ideal joint, discussed in Chapter 2, in that the pin is free to move within the inner perimeter of the bushing as dictated by the dyna mics of the system. Thus the centers of the components of the joint do not n ecessarily coincide. The ideal revolute joint was, however, modeled by assuming that the centers of th e joint components coincided at all times. Furthermore, the ideal revolute joint was simply considered to be a sing le point (axis) about which relative rotation between the two bodies occu rs. The actual physical joint was not modeled Two main approaches have been used to model the clearances in multibody systems. In one approach, the components of th e joint are assumed to be in c ontact at all times, that is, the pin and the bushing are in c ontinuous contact. Furuhashi et al [69] modeled the clearance as a fixed distance between the two centers of the joint compone nts. The fixed distance was interpreted as a massless link that introduced an extra degreeoffreedom. This technique has also been used by several other researchers [33, 38 ] to investigate the effect of joint clearance on the motion and dynamics of mechanical systems. In a different approach the components of th e revolute joint are not assumed to be in continuous contact. Rather, three di fferent configurations of the joint components are considered. In one of configurations, the pin is assumed to be in free flight [17]. This is when the pin is completely detached from the bushing and moves freely within the bushing depending on the system dynamics. In another configuration, the pin establishes contact with the bushing through an impact motion, and in the final mode the pi n moves while in contact with the bushing. The latter configuration has been re ferred to as following motion [17] This approach appears to be PAGE 46 46 the most realistic approach available. This latter approach will be used to model the nonideal joint and subsequently be integr ated with a wear pred iction procedure in order to determine the wear evolution in a joint. In this work the later approach will be termed as the general imperfect joint. General Imperfect Revolute Joint The general imperfect revolute joint closely models the actual joint in which the pin and bushing can assume three different configurations; i.e., freefall, impact or the following motion. In the real joint the motion of the pin is restrict ed within the bushing. Thus the pin moves freely within the inner perimeter of bushing until contac t is established. When contact/impact occurs, both bodies deform locally and a reaction force is developed in the contact region. The reaction force has the effect of restricting the pin within the bushing. This behavior is modeled in the general imperfect joint model by imposing a force constraint. A force, determined by a contact law, is applied every time contac t is established between the two joint components. The effect of the force is that the pin is prevented from leavin g the bushing inner perimeter. In what follows a discussion of the contactimpact force model will be presented thereafter followed by a discussion on modeling the general imperfect joint. ContactImpact Force Model A number of papers dealing with the subjec t including contact and impact of bodies in multibody systems have been authored. Two different approaches with regard to the treatment of this problem have been reported. These approaches are referred to as discontinuous and continuous. In the discontinuous approach, discussed by Khulief et al [12], the collision is assumed to be instantaneous. The dynamic analysis during this period is decomposed into two parts namely; pre and postcollision. The analysis is conducted as usual until impact/contact is detected. Momentum balance is then performed to determine the postcollision velocities. The PAGE 47 47 velocities in the analysis are then updated and the analysis is then resumed and performed as usual. Due to the momentum balance (and veloc ity updates), discontinuities are observed in the velocities. Although the method has been found to be efficient, the assumption of instantaneous impact becomes invalid if the duration of contact in the collision is large. This would limit the use of the method since the system configuratio n would have changed appreciably over that duration [43]. In the continuous approach, a c ontinuous force is used to describe the contactimpact force that results from the interaction of the two bodies. During the period of contact (including the impact), a force normal to the plane of contact is developed. The scenario is accounted for in the dynamic analysis by including the contact forces in the Differential Algebraic Equations of motion (as applied forces). The forces need to be updated at each time step, during the contact, to reflect the changing system configuration. In orde r to apply this method, it is necessary to know how the forces vary during the impactcontact peri od. Two main contact laws have been used in determining these forces in multibody systems, wh ich are the KelvinVoigt viscoelastic model applied by Kulief et al [44] and a contact force model with hysteresis damping applied by Lankarani et al [43], Ravn [13,15] and Flores et al [14,1618]. The KelvinVoigt viscoelastic model is a line ar contact force model in which a parallel springdamper assembly is used. The spring repr esents the elasticity of the joint components, whereas the damper accounts for the energy diss ipated in the process of collision. Energy dissipation is simulated by multiplying the rebound force by a coefficient of restitution. The model, as reported by Ravn [13], is expressed as, N rK loading F Ke unloading (31) PAGE 48 48 where K, and re are the elastic stiffness, deforma tion/penetration and coefficient of restitution, respectively. As an example, the contact force for this model is plotted in Figure 31(b). The force is plotted for the penetration shown in Figure 32(a) with a stiffness v alue of 6510 N/m. The plot of contact force verses the relative penetration is shown in Figure 32(c). It can be seen that during the rebound (unloading) the contact force is less than that in the com pression. The dissipated kinetic energy is also shown as the enclosed area. Although the KelvinVoigt model can provide co ntinuous contact force as well as accounts for energy dissipation during contact, it suffers fr om the difficulty of determining the elastic stiffness K. In addition, the linear rela tionship between the penetrati on and the contact force is a major simplification since the contact force w ould depend on other factors in addition to the material such as the geometry of the bodies in contact. The contact force model with hysteresis damping for multibody systems was first applied by Lankarani et al [43]. The model is an extension of the Hertz contact theory for spherical bodies in contact. It is based on the idea of dissipation of energy in terms of internal energy. According to Lankarani et al [43], the continuous contact force can be decomposed into an elastic and a damping term. The elastic term repr esents elasticity whereas the damping term is concerned with the dissipated energy during collision. The mo del can be expressed as, n NFKD (32) where K, D, and are the stiffness coefficient, damping coefficient the relative penetration and the relative penetration velocity. n is an exponent whose value is taken to be 1.5n Unlike the KelvinVoigt contact force model, the stiffness coefficient K is determined based on the geometry and material properties of the bodies in contact. The stiffness expression for colliding spheres i and j is given as, PAGE 49 49 1 24 3ij ij ijRR K RR hh (33) In Eq. (33) Ri and Rj are the radii of the spheres. The variables hi and hj are material properties which are dependent on the Youngs modulus E and Poissons ratio and can be computed as, 21 m m mhmij E (34) Hunt et al [45] proposed an expression for the damping coeffi cient written as, nD (35) where is referred to as a hysteresis damping f actor. Considering the kinetic energy of the system before and after the collision, Lankarani et al [43] determined an expression for the hysteresis damping factor. The expression which depends on the stiffness K, the initial penetration velocity () and a coefficient of restitution re is written as 231 4rKe (36) The contact force model with hyste resis can then be written as, 2 1.531 1 4r Ne FK (37) To demonstrate the model, the contact force is plotted in Figure 33(b) of the penetration shown in Figure 33(a) with a stiffness value of 6510 N/m1.5. The plot of contact force verses the relative penetration is shown in Figure 33c. Once again it can be seen that during the rebound phase the co ntact force is less th an that in the compression, indicating that energy has been lost. The dissipated kinetic energy is shown as the enclosed area. PAGE 50 50 It should be noted that the contac t force with hysteresis model (Eq. (37)) was derived for contact of spherical bodies. However, the concerne d in this work is th e contact in a planar revolute joint which corresponds to cylindrical joint com ponents. Contac tImpact models for cylindrical components have been proposed but these modes are quite difficult to use. For instance Dubowsky [1], proposed a contact force model for cyli ndrical bodies that can be expressed as, 3=l n1ij ij N nijijRRa F a FRR (38) where a is the length of th e shorter cylinder, Ri and Rj are the radii of th e cylinders, and other term are as previously defined. It is clear fr om this model that determining the contact force NF is complicated, requiring some iterative pr ocedures such as the NewtonRahpson method. Despite the difficulty of determining the contact force using Eq. (38), it has been shown [17] that there are no considerable g ains in using the cylindrical model instead of the spherical contact model (eq. (37)) with the cylinder radii used in pl ace of the s phere radii. Furthermore the cylindrical models do not account for dissipated energy. Thus in this work the contact model with hysteresis damping is used to determine th e continuous contact forc e in the revolute joint during collision of the components (i.e. pin and bushing). In addition to the contact force, the modeling of the imperfect joint can be enhanced by considering the friction as a resu lt of the contact between the jo int components. In this work, Coulomb friction model has been adopted. However, other models can also be used. The model, whose line of action is perpendicular to that of the contact force, is expressed as f kNFF (39) PAGE 51 51 In this expression NF is the contact force determined by Eq. (37) and k is the coefficient of kinetic friction and can be determined th rough experiments as described by Schmitz et al [46]. Modeling the General Imperfect Joint Model A discussion of how to determine the contac t forces between the jo int components of the general imperfect revolute joints was presented. In what follows, modeling of the imperfect joint will be presented. The discussion here follows closely the work of Ravn [13] and Flores [17]. Consider a revolute joint with clearance c (general imperfect joint) as shown Figure 31. The joint consist of two com ponent s namely; a pin and a bushing. For the sake of generality it is assumed that the pin and the bushing are rigidly attached to bodies i and j. The two bodies, referenced in the global c oordinates, are shown in Figure 34. Body coordinates xi yi and xj yj are fixed to the center of masses of the bodies A and E respectively. The bodyfixed coordinates are oriented at angles i and j relative to the global horizontal axis. When contact between the two bodies is established, a contact force (described by Eq. (37)) is generated for the entire contact period. F urthermore if there is relative motion between the joint components, a kineti c friction (described by Eq. (39)) is experience. In Figure 34, the point of contact C is defined as the center of contact region between the pin and the bushing. This point can be located using the eccentric vector e which is a vector connecting the bushing center D and the pin center B. At the time of contact the eccentric vector will point in the direction of the contact. This vector is expressed as, iiijjjerAsrAs (310) where ir and j r are vectors linking the global origin and the center or masses of the bodies. is and js are vectors in the local coordinate system that link the center of masses to the pin and PAGE 52 52 bushing centers respectively. iA and j A are transformation matrices that transform a vector from the local coordinate system to the global syst em. In this particular case they transformation vectors is and js into their global equivalent. The location of contact point C with respect to the pin and bushing can then be expressed as, ,C mmmmm R mijrrAsn. (311) In Eq. (311) Rm are the pin and bushing radii and n is a unit vector in the direction of e and is written as, e e e n e. (312) The penetration between the pin and bushing during the contact is computed as the different between the eccentric a nd the clearance as shown in Eq. (313). When the pin is not in contac t with the bushing the eccentric is smalle r than the clearance, and the penetration has a negative value. When the penetration has a va lue equal or greater than zero, contact is established. Thus when is greater than zero a contact force is applied between the bodies. The contact force vanishes when is equal to or less than zero. Thes e configurations are depicted in Figure 35. ec (313) In determining the contact force (using Eq. (37)) the relative penetr ation velocity is also requir ed. This will essentially be the difference between the velocities of the contact point. The velocity of each contact point is obtained as, ,mC mmmmmm R mijrr sn (314) PAGE 53 53 The relative velocities in the norma l and tangential direction can th en be computes as shown in (315). In Eq. (315), n is the unit tangent vector defined as nkn and k is the unit vector in the global zcoordinate. In Eq. (315), the normal velocity nv is also the relative penetration velocity and is positive during the penetration peri od and negative during the rebound period. CC nij CC tijv v rrn rrn (315) Once the relative penetration velocity has been computed, the contact and friction forces can be computed and assembled into the DAE to determine the dynamics of the system. It is noted that since the bodyfixed coordinates we re fixed at the center of masses of bodies i and j, the forces must also be applied at these locations rather than at the points of contact. Thus the transfer of the loads to the mass centers will result in an addition moment in each body. The force and moment equation for body ac ting at the center of mass for body i is given as, iN C iiii FFF mFrr, (316) where NN kNF F Fn Fn. (317) The corresponding loads for body j are jN C jjjj FFF mFrr. (318) Once the forces (contact and friction forces) for the nonidea l revolute joint are known, the description of the joint is co mplete. It should be noted that no kinematic constraint was PAGE 54 54 introduced while describing the non ideal joint. Instead a force constraint was used in the description. As a result a multibody system with th is type of nonideal revolute joint will have 2 additional degreesoffreedom whic h are catered for by the force constraint. Thus a kinematic analysis on a multibody system with a general imperfect joint (to determin e the motion of the system) is not possible. Instead the motion and the dynamics of the system must be determined through the dynamic analysis. A dem onstration of incorporating the general imperfect joint into a multibody system for a dynamic analysis will be presented in the following section. Example: SliderCrank Mechanism with Joint Clearance Between the Crank and the Connecting Rod (General Imperfect Joint) The use of the general imperfect joint model will be demonstrated by modeling a slidercrank mechanism that has a nonideal joint. Figure 36 shows a diagram of the sliderc rank mechanism which consists of four components (ground, crank, connecting rod and s lider). The components are connected to each other by three revolute joints a nd a translational joint. For this example, the revolute joint between the crank and the connectin g rod is modeled as a general imperfect joint. The dimension and mass properties for the mechanism are shown in Table 31. Also, the radii and material properties for the joint components (the pin and bushing) are shown in Table 32. For the contact force m odel, a value of 0.13 and 0.8 were used fo r the friction coefficient and the coefficient of restitution, respectively. These are summarized in Table 33. It is assumed in this analysis that the pin and the bushing are rigidl y attached to the crank and th e connecting rod, respectively. The crank is assumed to rotate at a constant angular ve locity of 30 rpm ( rad/sec). The kinematic constraint equations for this mechanism can be obtained using the procedures described in the previ ous Chapter. For this mechanism the constraint equations are as shown in Eq. (319). In Eq. (319) the first two rows describe the revolute joint between the PAGE 55 55 crank and the ground, whereas th e third and fourth describe the revolute join t between the connecting rod and the slider. The fifth and sixt h rows model the translational motion of the slider. The last row is the driving constraint that specifies the constant angular velocity of the crank. 111 111 2322 2322 3 3 1cos 0 sin 0 cos 0 sin 0 0 0 0 xl yl xxl yyl y t (319) It is noted that no constraint for the joint with clearance appears in Eq. (319). Instead a force constraint as previously desc ribed is used to restrict the m o tions of the crank relative to the connecting rod. The differential al gebraic equations of motion (Eq. (217)) repeated here as 22A T t q q qQ M q 0 q (320) can be assembled and solved (using techniques discussed in Chapter 2) to determined the dynamics of the system. A value of 250 is used for both stabilization parameters and Figure 37 shows representative results from the dynam ic analysis of the mechanism. In the figure a comparison of the jo int reaction force betw een the ideal and the nonideal joints is shown for four clearance si zes. In the first diagram th e clearance size is 0.5 m. For this case the two plots almost overlay each other as would be e xpected since the clearan ce is effectively zero. It should be mentioned that a value of zero for the clearance would lead to numerical difficulties since the eccentric vector will assume a value of zero (see Eq. (312)). As the clearance size is increased to 5mm the change in the system dyna mics becomes quite evident. The curve of the PAGE 56 56 reaction force is seen to evolve from a smooth one to one characterized by peak forces that increase in size with the clearance. In addition, th e general magnitude of the force is also seen to be higher for the case of the nonideal joint. The location of the peak contac t force can be explained by noti ng that they occur right after the minimum contact force (reaction force) be tween the joint components. This minimum corresponds to the slidercrank configuration when contact be tween the components is also minimal. Thereafter the contac t and thus the reaction force between the pin and the bushing increases rapidly. For the case of large clearance the pin has more freedom to move within the bushing and will therefore develop larger contact force peaks as the contact increases. It should be mentioned that the slidercrank configuration when there is minimal contact between the joint components is also a candidate fo r contact loss between the components. This is especially true when the clearance size is large. The last diagram in Figure 37 shows the c ontact force between the pin and the bushing when the clearance is quite la rge (5mm). The diag ram shows that contact between the two components is briefly lost when the slidercrank assumes the concerned configuration. Since the clearance is larger, the peak contact for ce that results after contact is reestablished is also higher than was previously observed. Summary and Discussion In this Chapter, a nonideal revolute joint m odel that closely represents a real joint was presented. The model, termed as the general im perfect joint, assumes th at the joint components (pin and bushing) locally deform /penetrate when they are in contact. The contact between the joint components generates a cont act force that is determined by a modified form of the Hertz contact law. This force prevents further penetra tion and restricts the motion of the pin within the inner perimeter of the bushing. Thus the imperfect joint is mode led using force constraints and not the kinematic constraint that we re used in the previous Chapter. PAGE 57 57 The developed joint model was used to model a nonideal join t in a slidercrank mechanism. In the mechanism, a clearance was pr esent in the revolute joint between the crank and the connecting rod. A dynamic analysis for this mechanism sh owed that the dynamics of the system is altered when the clearance is present at the joint. As the clearance is increased, force peaks are observed in addition to an increase in the force magnitude. These observations can be used to infer that wear at the joints also alters the system dynamics in a similar manner. PAGE 58 58 Table 31. Dimension and mass pa rameter for slidercrank mechanism Length (mm) Mass (g) Moment of inertia (kgm2) Crank 1.00 10.00 45.00 Connecting Rod 1.75 15.00 35.00 Slider 30.00 Table 32. Material proper ties for the joint components Pin Bushing Youngs modulus 0.29 0.38 Poisson ratio 206.8 GPa 0.5 GPa Radius 20 mm20.0005, 20.05, 20.5, 25 mm Table 33. Parameters for the force model Parameter Value Friction coefficient (pin & Bushing) 0.13 Coefficient of restitution (re) 0.8 PAGE 59 59 Figure 31. A revolute joint with clearance 0.331 0.433 0 0.1155 Time mm 0.331 0.433 0 600 TimeForce N 0 0.1155 600 mmForce N Figure 32. KelvinVoigt viscoelastic model. a) Relative pe netration; b) C ontact force during contact; c) Force verses the rela tive penetration during contact. 0.331 0.433 0 0.1155 Time mm 0.331 0.433 0 600 TimeForce N 0 0.1155 600 mmForce N Figure 33. Contactforce model with hystere sis damping (by Lankara ni & Nikravesh). a) Relative penetration; b) Contact force dur ing contact; c) Force verses the relative penetration during contact. Clearance c Bushing Pin (a) (c) (b) (a) (c) (b) PAGE 60 60 Figure 34. General imperfect joint. Figure 35. Penetrati on during contact between the pin and the bushing. 0 ec c0 ec 0 ec ePin Bushing e O y ri Si Sj rj A B E C D j i yi xi yj xj x e y1 x1 y2 x2 Bod y iBod y j PAGE 61 61 Figure 36. Slidercrank mech anisms with joint clearance. 8 10 12 14 16 18 0 1000 2000 3000 4000 5000 Reaction Force (clearance = 0.0005mm) Crack Angle (rad)Force (N) NonIdeal Joint Ideal Joint 8 10 12 14 16 18 0 1000 2000 3000 4000 5000 Crack Angle (rad)Force (N)Reaction Force (clearance = 0.05mm) NonIdeal Joint Ideal Joint 8 10 12 14 16 18 0 1000 2000 3000 4000 5000 Crack Angle (rad)Force (N)Reaction Force (clearance = 0.5mm) NonIdeal Joint Ideal Joint 8 10 12 14 16 18 0 1000 2000 3000 4000 5000 Crack Angle (rad)Force (N)Reaction Force (clearance = 5mm) NonIdeal Joint Ideal Joint Figure 37. Comparison of reac tion between the ideal and the nonideal joints for various joint clearance. Cran k Connectin g rod Slider Joint clearance Ground Ground y3 x3 x1 y1 y2 x2 y x PAGE 62 62 CHAPTER 4 WEARPREDICTION METHODOLOGY The procedu re to predict wear occurring at the interface of bodi es in contact and in relative motion will be presented in this Chapter. Emphasis is placed on ways to reduce computational costs while ensuring accuracy a nd stability in the predictions. Introduction When two bodies are in contact and are in relative motion, with respect to each other, wear is expected to develop on the regions of contact [47]. Similarly, in the case of a revolute joint wear is expected to occur since the two components of the joint (pin and bushing) are in contact and revolve relative to each other. Intuitively, it is safe to claim that the amount of wear at such a joint is affected by the type of material (which the joint components are made of), the relative sliding distance between the tw o components and the operating c onditions. Here, the operating conditions refers to the amount of reaction force developed at the joint and the condition of the joint which could be dry, lubricated, or contaminated with impurities. An enormous amount of effort has been placed into developing techniques to predict wear. A general trend that has emerged is the use of Archards wear model in an iterative procedure. Archards model is a linear wear model that estimates wear based on information about the contact condition (contact pressure and sliding di stance) and Tribological data that reflects information about the materials in contact and the operating conditions. Thus, wear procedures typically involve an iterative process in which in cremental wear is estimated at each iteration (based on the wear model) and accumulated up to the desired number of iterations. In earlier prediction procedures the linear model was employed to estimate worn geometry based on initial contact conditions The procedures assumed that the geometry and thus contact pressure did not evolve with wear. Thus only a single iteratio n was required in which linear PAGE 63 63 extrapolations were applied to determine the final geometry i.e. the geometry that would result after several thousand iterations. This proce dure has been found to produce erroneous results [48,49]. The primary reason for the inaccuracy is that the evolution of the geometry was neglected. In later procedures, wear pr edictions have been based on evolving contact conditions. The procedures allow the contact geometry to vary gradual and thus resulting in an iterative procedure in which the contact pressure and the sliding distance are comp uted at each iteration. The geometry is also updated at each iteration to reflect the worn material and thus geometry evolution. Depending on the effectiveness of the geometry update and the accuracy of the contact pressure calculat ions, the iterative procedures have be en found to yield useful results [5055]. This has allowed the procedure to be used in various application such as in the prediction of gear wear [5658], camfollower wear [5961], knee joint wear [62], and hi p joint wear [63,64]. Despite the reported success, the iterative procedure ha s been found to be quite computationally expensive. This is primarily due to the iterative process that is required to capture the evolving geometry. Several ideas have been implemented in an attempt to reduce computational costs associated with the wear simulation process. Pdra et al. [52] attempted to minimize the computational cost by using the Winkl er or elastic foundation model to determine the contact pressure distribution. The Winkler model was used as an alternative to the more expensive but relatively accurate FEM. Althoug h the method was found to be less expensive it can be argued that the benefit of using the more accurate results from the finite element technique outweigh the gains in computational efficiency when complicated geometries are considered. Pdra et al. [52] also employed a scaling appr oach to tackle the problem of computational costs. In this approach the incr emental wear at any particular cycle of the PAGE 64 64 simulation was scaled based on a predefined ma ximum allowable wear increment. The scaling factor was obtained as a rati o between the maximum allowable wear increment and the current maximum wear increment (maximum wear increment of entire ge ometry). They found that this procedure was more comput ationally effective. Kim et al. [50] used a constant extrapolation technique to reduce the computati onal costs for the wear problem. In their technique one finite element analysis was made to represent a number of wear cycles. Through this procedure, they were able to reduce the total number of analyses needed to estimate the final wear profile. A similar procedure was employed by McColl et al. [65] as well as Dickrell et al. [66]. In another paper [67], the computational costs of simula ting a pin on a rotating disc was reduced by approximating the state of strain on the center of the wear track as plain strain. A less costly twodimensional idealization was then used in place of the more expensive threedimensional problem. In what follows the widely used iterative we ar prediction procedure based on the Archards wear model will be presented. The procedure is then enhanced using an adaptive extrapolation scheme so as to reduce computational costs. The adaptive extrapolation scheme is an improvement to the constant extrapol ation scheme that was used by Kim et al. [50]. A pinpivot joint (a type of revolute joint), shown in Figure 41, will be used to demonstrate the procedure. W ear predictions on the pinpivot assembly will be conducted and validated with experiments. Wear Model As was mentioned the Archards wear law [ 68] forms the basis for the wearprediction methodology. In that model, first published by Holm [69], the worn out volume, during the process of wear, is considered to be proportional to the norm al load. The model is express mathematically as follows: PAGE 65 65 NF V K sH (41) where V is the volume lost, s the relative sliding distance betw een the two bodies in contact, K is the dimensionless wear coefficient, H is the Brinell hardness of the softer material, and NF is a normal force. Since the wear depth is the quantity of interest as opposed to the volume lost, Eq. (41) is usually written in the following form: NhA kF s (42) where h is the wear depth and A is the contact area such that the volume becomes VhA The nondimensioned wear coefficient K and the hardness are bundled up into a single dimensioned wear coefficient k (Pa1). It should be noted that the wear coefficient kis not an intrinsic material property but is also depend ent on the operating condition. Equation (42) can further be sim plified by noting that the contact pressure may be e xpressed with the relation N p FA so that the wear model is expressed as h kp s (43) The wear process is generally considered to be a dynamic process (rate of change of the wear depth with respect to sliding distance) so that the first order differential form of Eq. (43) can be expressed as: ()dh kps ds (44) where the sliding distance is considered as the time in the dynamic process, and the contact pressure is a function of the sliding distance. A numerical solution for the wear depth may be obtained by estimating the derivative in Eq. (44) with a finite divide differe nce to yield the depth as follows: PAGE 66 66 1jjjjhhkps (45) In Eq. (45), jh refers to the wear depth at the th j iteration while 1jh represents the wear depth at the previous iteration. The last term of Eq. (45) is the incremental wear depth which is a function of the contact p ressure and incremental sliding distance (js ) at the corresponding iteration. If information about the wear coefficientk, the contact pressure j p and the sliding distancejs is available at all iterations ( j ), the wear depth on a cont act interface fo r a specified sliding distance s can be estimated using Eq. (45). Here the sliding distance is an accumulation of the incremental sliding distance for all iterations (_niter) and is expressed as 1 niter j jss (46) The contact pressure ( p ) may be obtained using numeri cal methods. Various methods have been employed in computing the contact pr essure and sliding distance. These include the finite element method (FEM) [50,51,53,54,65,67,7072], boundary element method [73], Winkler model [52,56,74], Hertz contact model [57,58] each having its pros and cons with regard to accuracy and computational expe nse. In this work the FEM is used. The value of k for a specific operating condition and given pair of materials may be obtained by experiments [50, 75, 76]. It is worth noting that measured values of wear coefficients usually have large scatter and may affect wear predictions significantly. Care should thus be taken in obtaining these values. Uncertainty analysis for measured values of wear coefficients, such as those presented by Schmitz et al. [75], may be of considerable benefit. PAGE 67 67 Since all the necessary components of Eq. (45) can be determined, the simulations to predict wear can be conducted. A discussion on th e wear sim ulation procedure will be presented as follow. Wear Simulation Procedure A number of papers that demons trate the implementation of Eq. (45) in predicting wear, have been published [5052,54,56,65,77]. Although th e details of the various procedures differ, three m ain steps are common to all of them. These include the following: Computation of the contact pressure resulting from the contact of bodies. Determination of the incremental w ear amount based on the wear model. Update of geometry to reflect the wear am ount and to provide the new geometry for the next iteration. A brief discussion of these steps will follow. To facilitate the discussion, a pinpivot joint that experiences wear at the contact interface will be used. The pi npivot assembly is shown in Figure 41. It consists of a pin that oscillates inside a pivot. The objective woul d thus be to determ ine the amount of wear that occurs on both the pin and the pivot afte r several thousand pin oscillations. Computation of Contact Pressure The contact pressure, as previously mentioned, can be determined in a variety of ways. In this work, however, commercial finite element software, ANSYS, which is well adopted for this type of problem, is used. Figure 41 shows the diagram of the pinpivo t assem bly. The corresponding finite element model (a discretization of the pinpivot domain) is shown in Figure 42. Eightnode quadrilateral elem ents are used to model the pin and the pivot whereas threenode contact elements are used to represent the contact surface. The contact elemen ts are generated on the outer and inner surfaces PAGE 68 68 of the pin and pivot, respectively, and provide information about the contact between the two bodies. The use of the eightnode quadrilateral and threenode contact elements is not a limitation; other elements may also be used depending on the application and accuracy desired. Since the pinpivot model is discretized, the contact pressure on the contact surface (both pin and bushing contact surface) can be obtained at discrete locations, that is, at the nodal locations (or integration points). Details of how to perform the c ontact analysis can be found in the ANSYS users manual. Determining the Wear In addition to the spatial discre tization, the motion of the pin is also discretized in order to specify the incremental sliding distance (s ) (temporal discretization). A pin cycle which involves oscillation of the pin from one extreme to the other is disc retized into several steps. The discretization is such that each step corresponds to a specific pin angle between the two extremes from which the incremental sliding distance can be calculated. At each step a finite element analysis is perf ormed to determine the contact pressure over the contact region. The wear de pth during any cycle and at a ny location on the contact surface can then be determined by Eq. (47) which is a modification of Eq. (45) to include the discretizatio n of the cycles. ,,,1,, nijnijinihhkps (47) In Eq. (47), n refers to the contact surface node nu mbers which may or may not establish contact with the opposing surface. The subscript i and j indicate the current step and cycle, respectively. All other terms are as defined previously. The simulation procedure will involve dete rmining the incremental wear depth (, inikps ) at a particular step for the entire contact region. The geometry is then updated to reflect the wear PAGE 69 69 amount. The pin is then incrementally rotated to the next step and the corresponding incremental wear depth is computed. The geometry is once ag ain updated. This is repe ated until the cycle is completed. The next cycle iteration is started and the procedure repeated to the desired number of cycles. The simulation procedure is summarized in a flow chart shown in Figure 43. Geometry Update Procedure The process of geometry update is necessary in order to account for material loss through the process of wear. Indeed material removal changes the contact surface and causes a redistribution of the contact pressure. These change s can only be captured if the surface is altered through a geometry update. The geometry update pr ocedure involves two steps. These steps are outlined below: Determination of the normal direction (vector) of the contact surface at the location of each surface node (contact node). Shift the position of the surface nodes in the direction of the normal vector by an amount equal to the wear incr ement for each iteration. The contact surface is made up of a collection of contact elem ents. In the present case each contact element has three contact nod es. One such element is shown in Figure 44. Any location on the contact surface can be dete rm ined by an interpolating its coordinates from the surrounding nodal locations (surrounding contact nodes). Thus the contact surface can be expressed using nodal locations and a set of interpolation f unctions known as the shape functions. This information can be exploited to derive the ex pression for the normal di rection of the contact surface at each node. The shape functions for the contact element shown in Figure 44 may be written as: PAGE 70 70 1 2 31 1 2 11 1 1 2Ntt Ntt Ntt (48) where t is the local coordinate parameter. The surf ace of an element can then be described in terms of the nodal coordinates and as a function of the local coordinate. The expression for the surface is given in Eq. (49). 1 1 1232 1232 3 3000 000 x y NNNx x NNNy y x y (49) where k x and k y are the coordinates of node k (1,2,3k )for the element of interest. Thus a value of 1 t will supply the location of the first node 11(,) x y and a value of 1 t will supply the location of the third node 33(,) x y whereas the second node is determined by 0 t Thus if the nodal locations of the cont act nodes are known, any locati on of the contact surface is completely defined by expression in Eq. (49). If the vecto r tangent to the surface (c ontact element surface) is denoted as tv then its value for the element can be obtained as follows: 0xy tt tvi j k (410) where the partial differentials are given in Eq. (411) or (412). PAGE 71 71 1 1 3 12 2 2 3 12 3 3000 000 x y N NN x x tttt y yN NN t ttt x y (411) or 3 1 3 1 r r r r r rN x x tt N y y tt (412) The vector normal (nv ) to the surface can be expressed as a cross product of the tangent vector (tv ) and the vector perpendicular to the plane of the surface ( 0,0,1pv). This cross product is expressed in Eq. (413) where n denotes the contact node number. n n n tp n tpvv v vv (413) The resulting unit normal vector then appears as follows: 22 nyx tt x y tt ni j v (414) or, ,_,_, nnormxnnormynvv nvi j (415) where normxv and normyv are the components of the vector normal to the surface. Once the contact pressure dist ribution and normal vectors at all the nodes on the surface have been determined, the geometry can be updated. The update is done by moving the surface PAGE 72 72 nodes in the direction of the uni t normal vector. The coordinate of the new node position at any step of any cycle can be written as follows: ,, ,1, _, ,, ,1, _,nijnij normxn ni nijnij normxnxxv kps yyv (416) The process of the geometry update is shown in Figure 45. In this diagram the wear depth is gross ly exaggerated to illustrate the concept. The procedure for the geometry update has been used successfully in the wearsimulation process. Boundary Displacement The update procedure discussed in the previ ous section involved m oving only the nodes at the boundary to reflect the wear. If the am ount of w ear is large, the elemen ts associated with the boundary nodes will become distorted and the resu lts from the finite element analysis will become questionable. It is thus necessary to implement some strategy that will prevent this distortion. On possible solution is applying a procedure, termed as boundary displacement. The procedure, borrowed from structural shape optim ization [78], involves re positioning the internal nodes in addition to the boundary nodes. This ensures that the mesh regularity is maintained. Figure 46 illustrates how the internal node is reposition ed in order to reflect the wear. The process of repositioning the internal nodes in this approach can be broken down into three steps. The steps are outlined as follows: 1. Perform wear analysis and geometry update as has been discussed in the two previous sections. This is done until the amount of wear is equal to or exceeds a specified value. The wear analysis is then halted and the boundary displacement needs to be performed. 2. Using the initial geometry (befor e wear), perform a static fin ite element analysis using the accumulated wear depth as the displacement boundary condition. 3. Save the displacement of all the nodes. Usi ng this configuration of nodes as the updated geometry, proceed with the wear analysis (in step 1). PAGE 73 73 This procedure has previously be en used in wear analysis [50] and has been found prevent mesh distortion. In addition it allows for fine mesh to be used, allowing for more accurate finite element results. Reducing Computational Costs The procedu re discussed in the previous s ection provides a way to simulate the wear resulting from bodies in contact and in relativ e motion. However, the process can be quite expensive. For instance, if one desires to simu late 100,000 cycles for a case in which each cycle is discretized into 10 steps then 1,000,000 fi nite element analyses (nonlinear) as well as geometry updates would be required. Clearly th is may not be practical and the need for techniques to reduce the computational cost be comes immediately appare nt. One procedure that has previously [50,65,66] been employed to redu ce the computational costs is the use constant extrapolations. Extrapolations have been used in various forms with the goal of reducing computational costs. In this wo rk an extrapolation factor (A) is used to project the wear depth at a particular cycle to that of se veral hundreds of cycles Essentially, the extrapolation is the total number of cycles for which extrapolation is de sired. Thus according to this definition, the extrapolation factor can only take on positive integers values. Equation (47), repeated here as ,,,1,,nijnijinihhkps (417) can be modified slightly in order to incorporate an extrapolation f actor. It is first noted that the first term on the right hand side (R.H.S.) of Eq. (417) refers to the cumulative wear depth from previous cy cles whereas the last term refers to the incremental wear depth at the current step and cycle. As a way to minimize computationa l costs, it is assumed that the next A cycles (as many cycles as the value of the extrapolation) will have the same amount of wear depth as that of PAGE 74 74 the current step and cycle. The total incrementa l wear depth for those many cycles may then be obtained by multiplying the last term of Eq. (417) with the extrapolation factor. The resulting expression is shown in the following equation: ,,,1,,nijAnijinihhkAps (418) Utilizing the same concept, a new expression can be written to describe the position of the contact nodes for the geometry update proce ss. This expression is as follows: ,, ,1, _, ,, ,1, _,nijAnij normxn ni nijAnij normxnxxv kAps yyv (419) As may be expected, the level of accuracy of the wear simulation is reduced when extrapolations are used. This is directly rela ted to the assumption that the same value of incremental wear depth is maintained for se veral cycles. In realit y the geometry would continuously evolve resulting in a continuous redist ribution of the contact pressure. The evolving contact pressure would thus dictate a continuous change in the incremental wear depth at each cycle. Intuitively if the extrapolation sizes are small the corresponding inaccuracy will also be small. Use of extrapolations may also cause problems in simulation stability. Here stability is defined with regard to the cont act pressure distribution and hen ce the wear profile. An ideally stable wear simulation would be defined as one in which the contact pressure distribution remained smooth (with no sharp or sudden changes in the distributio n) for the entire duration of the simulation. It is however unlikely to have smooth pressure distribution throughout the simulation process. As a result a more relaxed definition of stability is adopted where by sudden changes in the pressure distribut ion are allowed to occur. In Figure 410a, the contact pressure in the contact interface of the pinp ivot assembly is seen to vary smoothly over the contact region PAGE 75 75 except for small peaks at the contact edges. The peaks are attributed to the transition from a region of contact to a region of no contact. This transition occurs at a point which can not be represented by a discrete model. The result is that there is an abrupt change in the surface curvature which causes high contact pressure. If su ch contact pressure dist ribution is maintained through out the simulation, the simulation can be referred to as stable. On the contrary, the diagram in Figure 410b is representative of contact pr es sure distribution that would constitute an unstable wear simulation. When very large extrapolation sizes ar e used, wavy pressure distributions ( Figure 410b) are obse rved and the simulation becomes unstable. The shift to instability due to the use of large extrapolation sizes can be explained as follows. The contact pressure dist ribution (obtained from the finite element analysis) is generally not pe rfectly smooth (e.g. the pressure peaks at the contact edges). This may be due to the discre tization error stemming form the finite element analysis. The use of an extrapolation factor ma gnifies these imperfections so that when the geometry is updated the contact surface smoothness is reduced. If large extrapolation sizes are used, the regions that experience high contact pres sure in a particular step of the simulation are worn out excessively so that in the following st ep these regions experience little or no contact. On the other hand, the regions that did not experience high contact pressure will be worn out less and thus will experience greater co ntact pressure in the next step This behavior will repeat in subsequent steps causing the surface to become increasingly unsmooth. The simulation will then become unstable. If, however, smaller extrapolati on sizes are used the wearing process acts as an optimizer to smoothen the surface but computational costs increase. This raises the question as to whether an optimum extrapolation size exists. In an effort to determine an optimum extrapolation PAGE 76 76 size, the adaptive extrapolation scheme was developed. This scheme will be discussed in the following section. Adaptive Extrapolation Scheme The adaptive extrapolation techni que is an idea proposed as an alternative to the constant extrapo lation scheme. The idea behind it is to se ek for the largest extrapolation size while maintaining a state of stability (smooth pressure distribution) th roughout the simulation process. The scheme is a threestep process. The first pa rt of the scheme involve s the selection of an initial extrapolation size (0A). The selection is based on e xperience, however it has been observed that as longs as the initi al value is not too la rge the wear prediction is not affected. This is because the adaptive scheme will automatically adjust to the appropriate value. In the second part of the adaptive extrapolati on scheme, a stability check is performed. The stability check involves monitoring the contact pressure distribution within an element for all elements on the contact surface. This essentially translates to monitoring the local pressure variation. If the contact pressure difference within an element is found to exceed a stated critical pressure difference crit p then a state of instability is noted and vi ce versa. In the final step of the adaptive scheme, the extrapolation size is altered based on the result of the stability check. That is, the extrapolation size is increased for the stable case and decreased for the unstable case. This process can be summarized as follows: 1incelecrit 1decelecrit if if j j jAApp A AApp (420) It must be mentioned that in order to mainta in consistency in the geometry update as well as in the bookkeeping of the number of simulated cycles, a single extrapolation size must be maintained through out a cycle. That is, every st ep in the discretized cycle will have the same PAGE 77 77 extrapolation size while different cycles may have different extrapolation sizes. Figure 48 shows a graph of the extrapolation history for the oscill ating pinpivot assem bly. From the graph, it can be seen that the extrapolati on took on a conservative initial value of about 3900 and increased steadily up to the 12th cycle (actual computer cycles not considering the extrapolations). Thereafter the extrapolation size os cillated about a mean of 6000. Experimental Validation Probably the m ost convincing wa y to validate the results of a simulation is to compare them against those from an actual experiment. In this work the simulation procedure is validated by comparing simulation results to results from a w ear tests. The wear test consisted of a fixed steel pin inside an unlubricated oscillating steel pivot as shown in Figure 41. The pivot w as set to oscillate with amplitude of 30 and was loaded in the direction of its shoulder. The resulting pressure at the crosssectional of the pivot was 60MPa. The load was kept approximately constant throughout the test. A total number of 408,000 cycles were completed during the test to yield a maximum wear depth of about 2mm. It should be noted, for the sake of comparison, that the definition of th e test cycles is differe nt from that of the simulation cycles. Here a test cycle is defined as a complete rotation from one extreme to the other and then back to the st arting position (in this case 30 to 30 and back to 30) where as in the simulation a cycle is considered to be a rota tion from one extreme to the other(i.e. case 30 to 30). The test information is summarized in Table 41. Sim ulation test were performed with the pin oscillation amplitude and pivot loading identical to that of the actual wear test The simulation test was run for 100,000 cycles (considering the extrapolation). A wear coefficient of 51.010 mm3/Nm (typical on unlubricated steel on steel contact) was used. This value is obtained from pinondisk tests results PAGE 78 78 reported by Kim et al. [50].The simulation test parameters are summarized in Table 42. In Figure 49, the history of wear for the pin and pivot nodes that experienced the m ost wear is shown. From the figure, a transient and steady stat e wear regime can be identified as discussed by Yang et al. [67]. The transient wear regime co rresponds to the beginning of the simulation. As the simulation progresses the geometry and the pressure distribution evolve and the wear transitions to a steady state regime. The evolut ion of the contact pre ssure is depicted in Error! Reference source not found. where contact pressures of nine in termediate cycles are plotted. Within the steady state wear regime, the wear is approximately linear with respect to the cycles as can be seen in Figure 49. This information may be exploited to approxim ate the wear after 408,000 cycles (since th e analysis was only done to 100,000 cycles). Through linear extrapolation, a value of 1.867mm was predic ted as the maximum wear depth on the pin. Although this value underestimates the wear depth it provides useful information about the wear. The variation of the extrapol ation size is depicted in Figure 411. Summary and Discussion In this Chapter a m ethodology to predict wear in a revolute joint was presented. The methodology is built upon a widely used iterative wear prediction procedure that is based on the linear Archards wear model. In the procedure, in cremental wear is estimated base on the contact pressure resulting from the contact of the joint components, the incremental sliding distance and the wear coefficient (which reflects both the material which the joint components are made of and the operating condition). The geometry is then updated so as to reflect the evolution of the surface and thus to account for changing contact cond itions. This process is then iterated and the incremental wear is accumulated up to the desired number cycles. To ensure that the computational costs are re duced, extrapolations are used. In the past extrapolations have been used to project the wear depth at a particular cy cle to that of several PAGE 79 79 hundreds of cycles and thus reducing computationa l costs. The difficulty associated with this approach is in determining an a ppropriate extrapolation size. This is because large extrapolation sizes would cause the geometry to distort after several iterations while small iterations would result in suboptimal use of resources. One of the contributions in this work was the development of the adaptive extrapolation scheme, which optim izes the selection of the extrapolation size. The adaptive scheme was developed to be used with the FEM. In the scheme, the extrapolation size is continuously increased until the geometry be gins to distort (or more specifically the mesh begins to distort). Distortion is monitored locally by comparing the pressure difference within an element, for all the surface elemen ts in the contact region, against a predefined critical. When the element pressure difference exceeds the predefined critical, distortion is noted and the extrapolation size is reduced. A validation of the wear procedure was al so conducted. The validation is done by comparing the results from the simulation to that of an experimental counterpart. The wear occurring at the contact interface of pinpivot assembly was simulated. The predicted wear depth deviated from the actual experimental wear depth by approximately 7%. Even though this deviation appears to be large th e predicted results is able to give a good insight into the wear occurring at the interface. Indeed like any other approximation tec hnique, errors are inherent. A number of factors contribute to this discrepancy including the wear model, which is not an exact representation of wear and the finite element analysis, which is an approximation technique. Another contributor is the wear coefficient. The wear coefficient is obtai ned experimentally and as was mentioned has a large scatter. Errors in th e wear coefficient considerably affect the results of the simulation. For instance, using a wear coefficient of 51.210 mm3/Nm resulted in predicted wear depth of 2.028mm. The new wear co efficient, which is still within the range of PAGE 80 80 scatter according to Kim et al. [50], has a deviation of about 1.4% from the experimental value. This is indeed a large improvement from the previ ous predictions. It is thus concluded that even though the procedure does not accurately predict th e wear the results obtain ed are of the correct order of magnitude and can be used for preliminary design. PAGE 81 81 Table 41. Wear test informa tion for the pinpivot assembly Test Parameters Values Oscillation amplitude 0 Load (crosssectional pressure) 60MPa Test condition Unlubricated steel on steel Total cycles 408,000 Max wear depth on pin ~2.00mm Table 42. Simulation parameters for the pinpivot simulation test Simulation Parameters Value Oscillation amplitude 0 Load (crosssectional pressure) 60MPa Wear coefficient (k) 51.010 mm3/Nm Total cycles 100,000 Steps per cycle 10 Table 43. Comparison of results form the simulation and Expt. wear tests Max. wear depth (pin) (mm) Simulation time (min.) Actual test 2.000 Step update 1.867 206 PAGE 82 82 Figure 41. Pinpivot assembly. Figure 42. Pinpivot finite element model. PAGE 83 83 Figure 43. Wear simulation flow ch art for the step update procedure. Figure 44. A threenode contact element used to represent the contact surface. In p ut Model Solve Contact problem (Contact Pressure) Application of Wear Rule 1) Determine wear amount 2) Determine new surface loc. U p date Model Total Cycles? End of Simulation C y cle Coun t Ste p coun t Total Steps/cycle ? PAGE 84 84 Figure 45. Geometry updates process. Figure 46. Geometry updates process using the boundary displacement approach. Figure 47. Contact pressure di stribution on a pinpivot assembly. a) The case of a stable wear simulation. b) The case of an unstable wear simulation. b a Initial Geometry Updated Geometry Old surface New surface PAGE 85 85 0 10 20 30 40 0 1000 2000 3000 4000 5000 6000 7000 CyclesExtrapolationExtrapolation Vs Cycles Figure 48. Extrapolation hist ory for a pinpivot assembly. 0 25 50 75 100 0 0.06 0.12 0.18 0.24 Cycles x 1000Wear Depth (mm)Wear on Pin & Pivot Pin Pivot Figure 49. Cumulative maxi mum wear on pin and pivot. PAGE 86 86 Figure 410. Evolution of the contact pr essure for nine intermediate cycles. PAGE 87 87 0 5 10 15 20 0 1000 2000 3000 4000 5000 6000 7000 CyclesExtrapolationExtrapolation Vs Cycles Figure 411. Extrapolation history plot for the step updating simulation procedure. PAGE 88 88 CHAPTER 5 INTEGRATED MODEL: SYSTEM DYN AMIC S AND WEAR PREDICTION The analysis of multibody systems with joint we ar will be presented in this Chapter. As before, the analysis will be re stricted to rigid multibody system s and only wear in the revolute joints will be considered. Introduction The procedure for analyzing m ultibody systems with joint clearance was presented in Chapter 3. The procedure enables information a bout the nonideal joint to be extracted. This information includes contact location, incremen tal sliding distance and joint reaction forces which are critical for wear prediction. Furthermor e the procedure is able to capture any changes in the system dynamics due to changes in the jo int dimensions. A wear prediction procedure to estimate the wear occurring in a revolute join t was also presented in Chapter 4. The wear prediction procedure can be integrated into the dynamic analysis of a multibody system in order to gain insight into the overall performance of a system when wear is pres ent at one or more of its revolute joints. The integrated model is comp osed of two parts namely; dynamic analysis and wear analysis. The model is discussed in the following subsections. Dynamic Analysis In the first part of the integr ation process a dynam ic analysis is performed to determine the joint reaction force and the incremental sliding distance. These are the two quantities required (from the dynamic analysis) to perform the wear analysis. The analysis is done for a complete cycle and the reaction force as well as the increm ental sliding distance is obtained at each time increment of the discretized range. The reaction force at the nonideal joint is determined by the contact force model. Th us for a nonideal joint b, the contact and friction force at time increment ti is expressed as: PAGE 89 89 ,,, ,,,iii iiib NtNtt b tkNttF F Fn Fn (51) where NF is the contact force and n is a unit normal vector pointing in the direction of contact. NF and n are describe in Eq. (37) and (312) respectively. The values of NF and n are computed during the integration of the equations of motion. It is worth noting that in the case of the ideal joints, the reaction force is determined by the product of the Jacobian and the Lagrange multiplier. For instance, the react ion force for an ideal joint k can be obtained as shown in the following expression: jbTb jrF (52) where b jF is the reaction force at the revolute joint k corresponding to body j, jT r and b are the Jacobian submatrix and the Lagrange multipliers corresponding to the concerned revolute joint. The incremental sliding distance is also obt ained at each increment and described as: 1 iiitjttsR (53) where j R is the bushing radius and it is the angle difference (in radians) between the local xaxes of the two bodies i and j that share a revolute joint at a current time whereas 1 it corresponds to the difference at a previous time. The value of t is obtained as ij Wear Analysis The second part of the integra tion process involves a wear anal ysis. The am ount of wear is determined at each increment based on the reacti on force and sliding distance from the previous dynamic analysis. The reaction forc e at each time increment is us ed to determine the contact pressure (through a finite element analysis) be tween the joint components (pin and bushing). Incremental wear can then be computed a nd the geometry is updated (see Chapter 4). PAGE 90 90 The wear prescribed by geomet ry update increases the cleara nce. The inner perimeter of the bushing is altered from its circular shape to one dictated by the wear. This implies that the value of c in Eq. (313) is no longer a constant value. Instead c depends on the location of contact C (as defined in Figure 34 and in Eq. (311)) and must also to be updated during the dynam ic analysis. The corresponding update for c is give by icR C bjrr (54) where i R is the original pin radius, and C j r and br are the position vector of the contact point C and the bushing center, respectively. The la st term on the left hand side of Eq. (54) ( C bjrr ) is the distance between the bushing center and the contact location C. This is the quantity that will vary as the bushing is worn out. After the pin and bushing geometry is updated an d the clearance size is adjusted to reflect the wear, another dynamic analysis is performed. The wear is th en computed based on the result of the new dynamic analysis. The geometry and the clearance size are once again updated. This process is iterated to the desired number of cy cles. The process is summarized in the flowchart shown in Figure 51. Demonstration of the Integration Process In this sec tion the integrated model will be used to predict the wear occurring at the joint of a slidercrank mechanism. The slidercrank mechan ism, encountered in Chapter 3, will be used to facilitate the demonstration. A diagram of the slidercrank is shown in Figure 52. The dim ensions and mass parameters for the slidercrank are shown in Table 51and Table 52. The dem onstration involves de termining the joint wear after several thousand revolutions of the crank. For this mechanism, wear is assume d to occur at the joint between the crank and the connecting rod. All other joints ar e assumed to have no wear. The jo int of interest consists of a PAGE 91 91 pin that is attached to the cr ank and a bushing attached to the connecting rod. To simplify the analysis, the pin is assumed to be made of ha rdened steel while the bushing is made of polytetrafluoroethylene (PTFE). This allows the pin to be modeled as a rigid body and the pin as a deformable body. In addition, the interaction of st eel and PTFE results in a very low wear rate for the steel and a high wear rate for PTFE. We ar on the pin is thus negligible and can be disregarded for low number of cycles. The con cern is therefore to acc ount for the wear on the PTFE bushing while performing the dynamic analysis. A value of 5.05x104 mm3/Nm, based on wear tests by Schm itz et al [75], was used for PTFE bushing wear rate. In additi on, the friction was considered at the joint and a value of 0.13 was used for the friction coefficient (Schmitz et al [46]). For the demonstration, the crank was operated at 30 rpm for 5,000 crank cycles. A summa ry of these simulation parameters is shown in Table 51 Representative results of the in itial dynamic analysis for the confi guration are shown in Figure 53 and Figure 54. A plot of the reaction for ce at th e joint of interest is shown Figure 53. The initial clearanc e for this joint was 0.005mm as listed in Table 52. In Figure 54 the locus of the cen ter of the contact region in the dynamic an alysis is shown. Thus each point on the graph corresponds to the center of the contact region. The center of the contact region was defined as the contact point C in Figure 34. From Figure 54b it can be seen that the loc ation of contact point C is concentrated on the left and the right side of the bushing inner diameter. This corresponds to bushing angular coordinates of approximately 0 and radians as defined in Figure 54c. This indicates that these two regions will expe rience more wear that other regions. PAGE 92 92 The wear result from the integrated mode l after 5,000 crank cycles is shown in Figure 55. The figure shows a plot of th e wear on the bushing as a func tion of the bushing angular coordinate. As was expected, two pe aks are observed to occur at 0 and radians. Thes e peaks correspond to the location where there was a concentration of the contact point C during the dynamic analysis (see Figure 54). The reason for the two p eaks can be explained by considering the m otion of the crank during a particular cycle. In the first half of the cycle, the crank tends to pull both the connectingrod and the slider wherea s in the second half of the cycle the crank pushes the rod and slider. This results in the pi n establishing contact wi th the bushing with the center of contact at approximately radians during the first half cycle and then establishing contact on the opposite side (0 radian s) during the second half cycle. Summary and Conclusions In this Chapter a procedure to analy ze m ultibody systems with revolute joint was presented. The procedure involves an alternatio n between dynamic analysis of the system and wear prediction of the concerned joint. The al ternating analyses are iterated to the desired number of cycles. Integration is achieved by updating the joint geometry based on the results of the wear analysis. The change s in the joint geometry are captured in successive dynamic analyses and thus the need for the alternation. The procedure is capable of determining the potential regions of contact and is able to quantify wear in thos e regions. The analysis of a slider crank mechanism with wear occurring at one of it s joints was presented to demonstrate the use of this procedure. PAGE 93 93 Table 51. Dimension and mass prope rties of the slidecrank mechanism Length (m) Mass (kg) Inertia x106 (kg.m2) Crank 1.00 10.00 45.00 Connecting rod 1.75 15.00 35.00 Slider 30.00 Table 52. Properties of the pin and bushing Pin Bushing Initial radius 20 mm 20.005m Depth 12.8mm Poisson ratio 0.29 0.38 Youngs Modulus 206.8 GPa 0.5 GPa Table 53. Test and simulation parameters Parameter Value Crank velocity 30 RPM ( rad/sec) Crank cycles 5,000 cycles Friction coefficient [46] 0.13 Wear coefficient [75] 5.05x104 mm3/Nm PAGE 94 94 Figure 51. Integration of wear anal ysis into system dynamics analysis Figure 52. Slidercrank mechanism with a w earing joint between the crank and the connecting rod. Cran k Connecting rod Slider Wearing Joint Ground Ground Input Geometr y Wear Anal y sis FEA Wear Rule Geometr y U p date D y namicAnal y sis Input Mechanism Parameters/Geometr y Assemble M, q AQ, Solve DAE for q and Integrate q for q and q Reaction Force Sliding Distance PAGE 95 95 7 8 9 10 11 12 0 1000 2000 3000 4000 5000 Force (N)Reaction Force (Contact Force Model) Crank Position (rad) Figure 53. Initial joint reaction fo rce for joint with clearance 0.0005mm. 30 15 0 15 30 30 15 0 15 30 Locus of the Center of Region Contact (Pin) xcoordinateycoordinate 30 15 0 15 30 30 15 0 15 30 Locus of the Center of Region Contact (Bushing) xcoordinate (mm)ycoordinate (mm) Figure 54. Locus of contact point C for a complete crank cycle. a) Locus for the pin. b) Locus for the bushing. c) Definition of the angular bushing coordinate. a b c x y PAGE 96 96 Figure 55. Wear on bushing as a function of th e bushing angular coordinate after 5,000 cycles. 4 2 0 2 4 0 0.1 0.2 0.3 0.4 0.5 Wear on Bushing Bushing Angular Coordinate (rad) ()Wear on Bushing (mm) FEM PAGE 97 97 CHAPTER 6 INTERGRATED MODEL: SYSTEM DYNAMIC S AND WEAR PR EDICTION USING THE ELASTIC FOUNDATION MODEL The use of the Elastic Foundation model in the analysis of planar multibody systems with joint wear will be presented in this Chapter. Introduction In Chapter 4 a procedure to predict wear occu rring at the interface of bodies in contact and relative motion was presented. The procedure requ ired determination of the contact pressure resulting from the contact of the bodies. This wa s achieved using the FEM (FEM). An alternative technique is the elastic foundation model (EFM) also known as the Winkler Surface model. This model has previously been used in wear pr ediction procedures. Flodin and Andersson [56] simulated mild wear in spur gears using th e EFM to compute contact pressure. Podra and Andersson [53] used the EFM to simulate the we ar for a sliding sphere on flat and cylinder on flat. They compared the wear results with result for wear simulations conducted with FEM. They reported reasonable agreements. In addition to computing the contact pressure the EFM may also be to compute the joint reaction force in a multibody dynamics framework. This approach was employed by Bei and Fregly [79] to perform multibody dyn amic simulation of knee contact. In this Chapter the use of the elastic foundation model in the analysis of multibody system will be presented. Elastic Foundation Model In the EFM the contact s urface is modeled as a set of springs (spread over the contact surface). The model is derived from plane strain elastic theory where an elastic layer of known thickness is bounded onto a rigid su rface [80]. The springs represen t the elastic layer and the thickness of the layer is composed of the thickness of one or both bodies (depending on whether one of the bodies is defined as rigid). The defo rmation on the elastic layer is produced by a rigid PAGE 98 98 indenter when a force Fin is applied on the indenter. The contact surface of the indenter takes into account the shape of the two bodies [80]. A figure of the EFM is shown in Figure 61. The EFM assum es that the springs are indepe ndent from each other and thus the shear force between them is neglected. The consequen ce of this assumption is that the model does not account for how pressure applied at one location affects the deforma tion at other location. This is contrary to what is experience in elastic c ontact where the displacement at one location is a function of the pressure applie d at other locations. Although this simplifying assumption violates the very nature of contact problems, some benefits can be derived from its use. In particular, the simplified model results in reduced pressure co mputational costs, faci litated analysis of conformal geometry, layered contact and nonlinear materials [79]. In addition, when the EFM is used in the analysis of multibody systems with jo int wear, the joint reactio n forces required for dynamic analysis can be recovered as a byproduct of the contact pressure. The need for a contact force model, therefore, disappears. This additional advantage will become clear when the discussion of the use of the EFM in multibody analysis is presented. The contact pressure for any spring (spring i) in the elastic foundation can be calculated from [56] W ii iE p L (61) where i p is the contact pressure, WE is the elastic modulus for the elastic layer, iL is the thickness of the elastic layer and i is the deformation of the spring. When both bodies are deformable WE is a composite of the elastic modulus and Poisson ratio of the two bodies. The procedure to determine the composite modulus is discussed by Podra [52] and in more detail by PAGE 99 99 Johnson [80]. For the purpose of illu stration, it is assumed that onl y one of the bodies in contact is deformable. For this case, a common expression for WE is give by [79, 8188] 1 112WE E (62) where E and are the elastic modulus and Poisson ratio of the deformable body, respectively. The contact pressure for the spring i can then be determined from 1 112i i iE p L (63) The total load supported by the elastic layer can then be computed as Load iiFpA (64) where iA is the area of each spring element. It shoul d be noted that total load supported by the elastic layer must equal the load applied on the indenter. Thus the followi ng equation must hold: InLoadFF (65) If the shape of the rigid indenter and the deformation of a pa rticular spring are known, the deformation at other springs can be determined. Equation (63) can then be used to determine contact pressure. It is, however, not typical th at the deform ation in springs is known. Rather, what is common is that the indenter force (shown as Fin in Figure 61) which causes the defor mation is known. The task is then to dete rmine the contact pressure that is caused by the force. This requires an iterative pr ocedure that is outlined as follows: 1. Guess an initial deformation for one of spring. 2. Based on the shape of the indent er and the initial guess, determine the deformation for the other springs. 3. Compute the contact pressure using Eq. (63). 4. Find the sum of the load (LoadF ) supported by the elas tic layer using Eq. (64). PAGE 100 100 5. Check if Eq. (65) is satisfied within some specified tole rance. If it is, stop the iteration otherwise repeat part 14 until convergence is achieved. The iteration procedure is summarized in Figure 62. Analysis of Multibody Systems with Joint Wear Using the EFM In this section, a discussi on of how the EFM can be used in the analysis of multibody systems with imperfect revolute joints will be pr esented. While the framework for this analysis, presented in Chapter 5, remains the same, the manner in which the revolute joint is modeled will be different. In addition the w ear prediction procedure will in volve the EFM rather than the FEM. Modeling the NonIdeal Revolute Joint Using the EFM The imperfect revolute joint was defined as consisting of two components, a pin and a bushing, that are rigidl y attached to bodies i and j which share the joint. A realization of the revolute joint with the two bodi es, referenced in the global coordinates, is shown in Figure 31. The im perfect revolute joint was modeled so that the pin and the bushing centers do not necessarily coincide during the mo tion of the mechanism. The pin is allowed to move within the bushing inner perimeter. Whenever contact is established, a contact force, determined by Eq. (37), is applied so as to rest rict the pin within the bushing. Instead of us ing the contact force model, the EFM can be used to de termine the appropriate force required to restrict the pin within the bus hing inner perimeter. The procedure will first require determination of the c ontact pressure resulting from the contact between the pin and bushing using Eq. (63). In this equation, the material properties ( E and ) and the thickness of the elastic layer (iL ) will generally be known from the pr oblem definition. The deformation, however, (i ) is dependent on the system dynamics. Referring to Figure 63, the deformation m ay be computed as: PAGE 101 101 ec (66) where e is the eccentricity and is defined as the distance between the pin and the bushing, and c is the initial clearance in th e joint. When the pin is not in contact with the bushing the eccentricity is smaller than the clearance, and th e penetration has a nega tive value. When the penetration has a value equal or greater than zero, contact is esta blished. Thus when is greater than zero a contact pressure can be computed using Eq.(63). The contact pressure vanishes when is equal to or less than zero. Th ese configurations are depicted in Figure 64. The eccentricity can be d etermined from the eccent ric vector which is a vector that points in the direction potential contact. As previously discussed in Chapter 3, the vector is calculated from iiijjjerAsrAs, (67) where ir and j r are vectors linking the global origin and the center of ma sses of the bodies. is and js are vectors in the local coordinate system that link the center of masses to the pin and bushing centers respectively. iA and j A are transformation matrices that transform a vector from the local coordinate system to the globa l system. The eccentricity is then give as the magnitude of the eccentric vector defined as e e. (68) Once the contact pressure has been evaluated, the contact force can then be computed from Eq. (64). It is also possible to compute the fricti on force using the Coulom b friction model discussed in Chapter 3. Both the contact and the friction forces may then be assembled into the differential algebraic equations of motion whic h can then be solved (as discu ssed in Chapter 2) to determine the system dynamics. It should be noted that si nce the bodyfixed coordinates were fixed at the PAGE 102 102 center of masses of bodies i and j the forces must also be applied at these locations rather than at the points of contact. Thus the transfer of the loads to the mass centers will result in an additional moment in each body. The appropriate tr ansformations are expressed in Eqs. (316), (317) and (318). Example: dynamic analysis of a slidercra nk with an imperfect joint using EFM. The slidercrank mechanism encountered in Chapters 2 and 3 will be used to demonstrate the use of the EFM in analysis of multibody system with im perfect joints. A diagram of the slidercrank mechanism is shown in Figure 65. The revolute joint betwee n the crank and the connecting rod is m odeled as an imperfect joint. All other join ts, including the translational joint (slider), are modeled as ideal. The dimension and mass pr operties for the mechanism are shown in Table 61. For this ex ample, it is assumed that the pin is made of steel, while the bushing is made of PTFE. The material properties as well as the ra dii of the pin and bushing are listed in Table 62. The crank is constrained to ro tate at a constant angul ar velocity of 30 rpm ( rad/sec). The kinematic constraint equations for this mechanism are identical to the one discussed in Chapter 4 and are repeated in Eq. (69) for convenience. It is worth noting that only driving constraint and the constrai n t equations representing the ideal joints appear in the set of kinematic constraint equations. The dynamics of the slider crank mechanism can then be determined by assembling and solving the differential algebraic equations of motion. 111 111 2322 2322 3 3 1cos 0 sin 0 cos 0 sin 0 0 0 0 xl yl xxl yyl y t (69) PAGE 103 103 Representative results of the dynamic analysis are shown in Figure 66, Figure 67 and Figure 68. In Figure 66 the joint reaction force (between the crank and th e connecting rod) for the EFM based model is compared to the joint fo rce from the ideal joint model. For the nonideal model a clearance of 0.0 005 mm was used. Since this clearance is quite small, it is expected that the response of the EFM based model should closel y resemble those from the ideal model as is observed in Figure 66. In Figure 67, the locus of the center of the contact region is plotted for both the pin and the bushing. From this plot it is seen that th e pin experiences contact on only one side whereas the bushing has a concentratio n of the contact points at two locations. This behavior was observed in the cas e of the previous analysis pr ocedure based on the Hertzian contact model (see Figure 54). Comparison of the predicte d f orce between the ideal model and the EFM base model for various joint clearances is shown in Figure 68. It can be seen from the plots th at as the clearance incr eases, the reaction force begins to exhibit some oscillations. The amplitude and frequency predicted by the EFM m odel differ from those predicted by contact force model (see Figure 37). However the trend, in term s of t he location of the oscillation is quite similar. Integrated Model: System Dynamics and Wear Prediction Using EFM It was seen, in the previous section, that in order to perform dynamic analysis of a system with imperfect revolute joints using the EFM, it was first necessary to determine the contact pressure at the revolute joint. Contact pressure information is thus available after every dynamic analysis and can be used for wear prediction using the procedures discussed in Chapter 4. The integration of the wear predic tion procedure into the dynamic an alysis is achieved by updating the bushing geometry after every wear analysis1. Since the inner perime ter of the bushing is 1 For simplicity it has been assumed that one of the joint components (pin) is rigid and thus does not wear. PAGE 104 104 altered from its circular shape to one dictated by the wear, the value of c in Eq. (66) is no longer a constant value but instead depends on the location of contact C (as defined in Figure 63 and in Eq. (311)). Thus, the value of c must also be updated in progressive dynamic analyses. The expression used to update the clearance c was discussed in Chapter 3 and is repeated here for convenience as: icR C bjrr (610) where i R is the original pin radius, and C j r and br are the position vector of the contact point C and the bushing center, respectively. The first term on the left hand side of Eq. (610) ( C b j rr) is the distance between the bushing center and the contact location C This is the quantity that will vary as the bushing is worn out. Example: analysis of a slidercrank with joint wear using EFM. Wear prediction at a joint of the slidercrank mechanism will be used to demonstrated integrated model based on EFM. For purposes of comparison with integrat ed model based on FEM, simulation parameters are selected similar to the example discussed in Chapter 5. The slidercrank mechanism is shown in Figure 65 and its dimensions a nd m ass properties are listed in Table 61. The pin and the bushing are assum ed to be made of steel and PTFE respectively. The material properties of the pin and bushing are listed in Table 62 and their radii and bushing depth are listed in Table 63. The crank is constrained to revolve at a constant angular of 30 rpm for 5,000 cycles. A plot of the locus of the center of the contact region (contact point C as defined in Figure 63) for the initial dynamics is shown in Figure 67. It can be s een that for th e pin ( Figure 67a), the contact point C is concen trated on the left side whereas for the bushing ( Figure 67b) the contact p oint C is concentrated on both the left a nd the right side of the bushing inner diameter. This corresponds to bushing a ngular coordinates of approximately 0 and radians as PAGE 105 105 defined in Figure 67c. This indicates that these tw o regions will expe rience more wear than other regions. The same behavior was observed in the integrated model based on the FEM (see Figure 54). In Figure 69 and Table 64, the wear predicted by the integrated model base on the EFM after 5,000 crank cycles is com pared to the wear predicted by the FEM based model. In the figure, wear is plotted as a f unction of the bushing angular co ordinate. As was expected, two peaks are observed to occur at approximately 0 and radians. The maximum wear depth predicted by the FEM based model is seen to be greater than that of FEM based model. On the contrary, the EFM based model predicted more wear in the region where the FEM based model experienced the least amount wear. This differenc e in wear profile is directly linked to the assumption in the EFM that the individual spri ngs are independent of each other and lateral effects are neglected. Another interesting observati on from wear results, in Figure 69 and Table 64, is that while the m aximum wear depth predicted by the two models differ, the wear profile for the models is such that the worn volumes predicted by the two models are approximately equal. This equality is a manifestation of the equality in the force as shown in Figure 66. It should be m ention however that as the wear increases in th e joint the force predictions for the two models differ slightly. As a result the wear volume will also differ. Also in Table 64 the computation time for the tw o m odels is compared. The EFM has a shorter computational time. This is because of two reasons: 1) the EFM is quite simple because the springs in the elastic layer are assumed to be independent, and 2) in the case of EFM based model only one set of analyses for each cycle is required to simultaneously determine both the contact force and the contact pressure during the dynamic analysis whereas in the case of the PAGE 106 106 FEM based model the contact force is determine in the dynamic analysis which in turn is used to determine the contact pressure. Summary and Conclusion In this Chapter the procedure to determine c ontact pressure and joint reaction force of a multibody system using the elastic foundation mode l was discussed. Wear prediction was then introduced into the dynamic analysis to allow for the analysis of the system with joint wear. A slidercrank mechanism was used to demonstrate the integrated model as well as for comparison with the FEM based model. The force prediction of the two models was id entical for near zero clearance but differed in magnitude as the clearan ce was increased. This difference is attributed to the different techniques employed in the models used to compute the force. It was seen that the maximum wear depth and the wear profile pr edicted by the two models differed. However, the wear volume predictions from the two mode ls were similar. The computation time for the EFM based integration model was found to be sh orter than that of the FEM based model. While it is well known that the FEM is more s uperior to the EFM, th e use of the EFM may be justified because of comput ational speed. The EFM based integrated model may be used to provides some qualitative information about the system. However results from the model should be interpreted with caution to avoid erroneous conclusions. PAGE 107 107 Table 61. Dimension and mass pa rameter for slidercrank mechanism Length (mm) Mass (g) Moment of inertia (kgm2) Crank 1.00 10.00 45.00 Connecting Rod 1.75 15.00 35.00 slider 30.00 Table 62. Material proper ties for the joint components Pin Bushing Youngs modulus 0.29 0.38 Poisson ratio 206.8 GPa 0.5 GPa Radius 20 mm20.00, 20.05, 20.5, 25 mm Depth 12.8 mm Friction coefficient (pin & Bushing) [46] 0.13 Table 63. Wear simulation parameters Parameter Value Crank speed 30 rpm Crank cycles 5,000 cycles Bushing inner diameter 20.005 mm Pin diameter 20.000 mm Bushing depth 12.8 mm Wear coefficient [75] 5.05x104 mm3/Nm Table 64. Comparison of results from prediction based on FEM and EFM FEM EFM Error Maximum wear depth 0.34 mm 0.29mm 13.0% Wear Volume 231.42 mm3241.76 mm3 4.3% Computation Time 162 min92 min PAGE 108 108 Figure 61. Elastic foundation model. Figure 62. Procedure to determine the contact pressure using the EFM. L Elastic Layer Rigid Surface i F I n Rigid Indenter Spring i End of Simulation Initial guess for 1 Determine i based shape of indenter and 1 Calculate Contact Pressure pi Find Fload L oad ii F pA Update value of 1 Fload= Fin ? PAGE 109 109 Figure 63. General imperfect joint. Figure 64. Penetrati on during contact between the pin and the bushing. Figure 65. Slidercrank mech anisms with joint clearance. 0 ec c0 ec 0 ec ePin Bushing e Cran k Connectin g rod Slider Joint clearance Ground Ground y3 x3 x1 y1 y2 x2 y x O y ri Si Sj rj A B E C D j i yi xi yj xj x e y1 x1 y2 x2 Bod y i Bod y j PAGE 110 110 8 10 12 14 16 18 0 1000 2000 3000 4000 5000 Crack Angle (rad)Force (N)Reaction Force (clearance = 0.0005mm) Elastic Foundation Model Contact Model Figure 66. Comparison of r eaction force between the ideal joint model and the EFM. 30 20 10 0 10 20 30 30 20 10 0 10 20 30 Locus of the Center of Region Contact (Pin) xcoordinateycoordinate 30 20 10 0 10 20 30 30 20 10 0 10 20 30 Locus of the Center of Region Contact (Bushing) xcoordinate (mm)ycoordinate (mm) Figure 67. Locus of contact point C for a complete crank cycle. a) Locus for the pin. b) Locus for the bushing. c) Definition of the angular bushing coordinate. a b c x y PAGE 111 111 8 10 12 14 16 18 0 1000 2000 3000 4000 5000 Crack Angle (rad)Force (N)Reaction Force (clearance = 0.0005mm) Elastic Foundation Model Ideal Model 8 10 12 14 16 18 0 1000 2000 3000 4000 5000 Crack Angle (rad)Force (N)Reaction Force (clearance = 0.05mm) Elastic Foundation Model Ideal Model 8 10 12 14 16 18 0 1000 2000 3000 4000 5000 Crack Angle (rad)Force (N)Reaction Force (clearance = 0.5mm) Elastic Foundation Model Ideal Model 8 10 12 14 16 18 0 1000 2000 3000 4000 5000 Crack Angle (rad)Force (N)Reaction Force (clearance = 5mm) Elastic Foundation Model Ideal Model Figure 68. Comparison of reaction force between ideal joint model EFM. PAGE 112 112 4 2 0 2 4 0 0.1 0.2 0.3 0.4 0.5 Comparison of Wear on Bushing Using FEM and EFM Bushing Angular Coordinate (rad) ()Wear on Bushing (mm) Finite Element Method Elastic Foundation Model Figure 69. Comparison of joint wear be tween the EFM and FEM after 5,000 cycles. PAGE 113 113 CHAPTER 7 EXPERIMENTAL VALIDATION OF THE INTEGRATED MODELS Introduction Two procedures to analyze multibody systems with joint wear were presented in Chapters 5 and 6. In addition predictions from the two pr ocedures were compared in Chapter 6. Before these models can be utilized for practical purpos es, their validity should be assessed. In this Chapter, the integrated models w ill be validated through experiments. Experiments for Model Validation For the purpose of validation, an experime ntal slidercrank mechanism was built. The mechanism was built at the Tribology Laboratory in the University of Florida. Likewise, all the wear tests were conducted at this laboratory. A diagram of the slidercrank mechanism, built for this purpose, is shown in Figure 71. The focus of the experimental tes t is to determ ine the wear that occurs at the joint between the crank and the connecting rod after several thousand revolutions of the crank. This joint essentially consists of a pin that is attached to the crank and a bushing attached to the connecting rod. The pin is made of hardened steel and is assumed to be hard enough so that no appreciable wear occurs on its surface. The bushing on the othe r hand is made of polytetrafluoroethylene (PTFE) which is soft and will experience consider able wear. A spring is attached to the slider which serves as a means to increase the joint reaction force and hence accelerates the wear occurring at the joint. The spring also ensures that contact between the pin and the bushing is maintained even in the presence of excessive wear. The slidercrank mechanism was built so as to minimize friction and wear (to a negligible amount ) at all joints except at the joint of interest (joint between crank and conn ecting rod). This is achieved by building the joint between connecting rod and the slider with a thrust air bearing and using a doveta il air bearing slider. PAGE 114 114 Details of the construction of the mechanism can be found in the works of Mauntler et al [89, 90]. The dimensions and mass properties for th e experimental slidercrank are shown in Table 71. The dim ensions, the material properties, the fi ction coefficient ant the wear coefficient of the joint components (bushing and pin) are listed in Table 72. The validation involves com paring the wear on the bushing from the experiment to the wear prediction by the two models (model based on FEM and EFM). For the validation, both the maximum wear depth and the wear profile (location of wear) will be compared for the two models. The comparison will thus give an indi cation of the performance of each model with respect to each other and with respect to the experiment. In the experimental test, the slidercrank mechanism was operated for 21,400 cycles. The crank was constrained to revolve at a constant velocity of 30 rpm. A spring with spring constant of 525 N/m was used. The test and simu lation parameters are summarized in Table 73. Figure 72 and Figure 73 show representative results from the initial system dynamics. Three plots of the joint reaction force from 1) the experiment 2) the model based on the EFM and 3) the model based on the contact force model are shown in Figure 72. The two models, which are id entical in this case, pred ict the joint reaction force reasonably well over the entire crank cycle except at about radians. At this location, the m easured force exhibits high frequency oscillation for a short duration. The location of these oscillations correspo nds to one half of the crank rotation when the slider changes direction. It is the belief of the author that these higher order dynamics is a result of the change in the di rection of the slider which most likely involves a slight rotation of the slider and thus a moment of brie f impact with the sliding rail. It should be mention that, although the magnitude of these osc illations is large, th ere effects on the wear PAGE 115 115 prediction is quite small. This is because the corres ponding incremental sliding distance is also quite small. Thus according to equation (45), repeated here for convenience as 1jjjjhhkps (71) the value of pj will be larger because of the large amplitude oscillations but the value of js will be quite small since the oscillation occur over a s hort duration. The result is that the incremental wear depth, for this region, will not vary signif icantly from the incremental wear depth of its neighboring regions. In Figure 73 the locus of the center of the re gion of contact are plotted for the pin and bushing. Figure 73(a) shows the locus of points when the model based in the EFM was used. It is seen that the entire pin surface will experience contact. On the othe r had the center of the region of contact on the bushing is concentrated on the left side of the bushing. This means that only one side of the bushing will experience wear and that the maximum wear will occur where there is a concentration of the center of the region of contact. The concentration of the contact point in this location is reasonable because the spring that restrict the motion of the pin relative to the bushing. Figure 73(b) corresponds to the plot of the locus or the center of the region of contact when the contact force m odel is used. It is clear that the two models give approximately identical predictions. The wear predicted by the FEM and the EFM models are compare in Figure 74 and Table 74. In Figure 74, it can be seen that while the FE M base m odel predicts a larger maximum wear depth, the EFM has a wider base. The wider ba se means that a wider region in the bushing surface is worn out. An interesting observation is that while the wear depth for the two models differs, their wear profile is such that the worn volume is equal for the two models. This equality PAGE 116 116 is a manifestation of the equality in the force as seen was in Figure 72. This behavior was identified an d explained in the previous chapter. The computation time for the wear predicti on based on the FEM and the EFM models is compared in Table 74. As was seen in the previous chapter the EFM model has a faster com putation time. The wear results from the experiment and the simulation results from the two models are compared in Figure 75, Table 75 and Table 76. From Figure 75(a) and Table 75, it can be seen that the m aximum wear dept h, the wear profile and the wear volume from the experiment have accurately been predicted by the FEA based model. There is however a discrepancy between 4.5< <6.3 that is attributed to the meas urement of the wear on the bushing in the experiment. In the case of the EFM model, the wear profile and the maximum wear depth are incorrectly predicted as shown in Figure 75(b) and Table 76. The location of maximum wear and the wear volum e are however correctly predicted as expected. Summary and Conclusions The objective of this Chapter was to validat e the procedures for analysis of multibody systems with joint wear. The two procedures we re presented in Chapters 5 and 6. For the validation, results from an experimental slidercrank mechanism were compared to results from the two models. In addition the performance of the two models was analyzed by comparing results between models and against the experimental results. The focus of the comparison was the joint force, contact location and wear on th e bushing at the joint be tween the crank and the connectingrod. All other joints were assumed to be free of wear and free friction as necessary; provisions were made for these assumptions. PAGE 117 117 For the initial dynamics the two models prov ided reasonably accurate prediction of the contact force. The two models pr oduced identical prediction fo r the location of maximum wear. Through the experiment, this location was later ve rified to be correct. Although this location was correctly predicted by both models, only the FEM based model gave an accurate prediction of the wear profile and maximum wear depth. Pr ediction from the FEM based model differed by 6.7% from the experiment, whereas the predic tion from the EFM model differed by 12.1% from the experiment. It should, however, be noted while the predictions on the wear profile differed, the wear volume predictions of the two models were identical and reasonably close to the experimental wear volume (8.2 %). Despite the poor prediction on the wear profil e and maximum wear depth, the EFM based model had a shorter computation time. The experiment took about 12 hours (excluding construction and setup time) and FEM based model took 11 hours while the EFM based model took 7 hours to complete. The speed of the EFM ba sed model is associated with the assumption that the springs in the elastic foundation model are independent. From the comparison between the results of th e two models and between the results of the two models and the experiments the following conc lusions can be made: 1) the two procedures (FEM and EFM based procedures) can accurately predict the contact force, th e contact locations and the wear volume, 2) FEM based procedure is a better predictor of maximum wear depth and the wear profile than the EFM based procedure, and 3) EFM model is computationally faster than the FEM based model. It can be concluded that the FEM based procedur e will be a better procedure in the analysis of multibody systems with joint wear when the co mputational costs in not of concern. On the other hand if the cost is of concern and only qua litative information about the system is needed PAGE 118 118 then the EFM based procedure will be the su itable choice. Other scenarios will require a compromise on either accuracy or computational costs. PAGE 119 119 Table 71. Dimension and mass prope rties of the slidecrank mechanism Length (m) Mass (kg) Inertia x106 (kg.m2) Crank 0.0381 0.4045 204.0 Connecting rod 0.1016 0.8175 5500.0 Slider 8.5000 Table 72. Properties of the pin and bushing Pin Bushing Bushing Inner radius 9.533 mm Outer radius 9.500 15.875mm Depth 13.100mm Poisson ratio 0.29 0.38 Density 7.8 g/cm32.25g/cm3Youngs Modulus 206.8 MPa 0.500 MPa Friction coefficient (steel & PTFE) [46] 0.13 Wear coefficient (steel & PTFE) [75] 5.05x104 mm3/Nm Table 73. Test and simulation parameters Properties Value Crank velocity 30 RPM Crank cycles 21400 Spring constant 525.4 N/m Table 74. Comparison of wear results fo r FEM and EFM models (21,400 crank cycles) Table 75. Comparison of wear results betw een test and FEM model (21,400 crank cycles) FEM EFM Difference Wear Volume 106.71 mm3106.68 mm3 0.02% Max wear depth 0.4779 mm0.4263 mm 10.70% Computation time 11hrs 7hrs 4 hrs Experimental Simulation (FEM) Difference Worn mass 0.2616 g 0.2401 g 8.2% Wear Volume 116.27 mm3106.71 mm3 8.2% Max wear depth 0.48500.4779 mm 1.5% PAGE 120 120 Table 76. Comparison of wear results betw een test and EFM model (21,400 crank cycles) Experimental Simulation (EFM) Difference Worn mass 0.2616 g 0.2400 g 8.2% Wear Volume 116.27 mm3106.68 mm3 8.2% Max wear depth 0.48500.4263 mm 12.1% PAGE 121 121 Figure 71. Experimental slider crank mechanism 0 1 2 3 4 5 6 0 50 100 150 200 Crank Position (rad)Joint Force (N)Comparison Joint Force Between Models and Expt EFM Contact Model Experiment Figure 72. Comparison of the initial joint reaction force between the two models and the experiment. PAGE 122 122 Figure 73. Locus of the center of the contact region. a) Pred iction based the elastic foundation model. b) Prediction based on the contact force model. 0 1 2 3 4 5 6 0 0.1 0.2 0.3 0.4 0.5 Comparison of Wear on Bushing Using FEM and EFM Bushing Angular Coordinate (rad) ()Wear on Bushing (mm) FEM EFM Figure 74. Comparison of the w ear prediction between the models. a 10 5 0 5 10 10 5 0 5 10 Locus of the Center of Region Contact (EFM) xcoordinate (mm)ycoordinate (mm) Pin Bushing b 10 5 0 5 10 10 5 0 5 10 Center of Region Contact (Contact Force Model) xcoordinate (mm)ycoordinate (mm) Pin Bushing x y PAGE 123 123 Figure 75. Comparison of the wear profile for the models a nd the experiment. a) Comparison between experiment and FEM. b) Comparison between experiment and EFM. a 0 1 2 3 4 5 6 0 0.1 0.2 0.3 0.4 0.5 Comparison of Wear (Expt and FEM) Bushing Angular Coordinate (rad) ()Wear on Bushing (mm) Expt FEM b 0 1 2 3 4 5 6 0 0.1 0.2 0.3 0.4 0.5 Comparison of Wear (Expt and EFM) Bushing Angular Coordinate (rad) ()Wear on Bushing (mm) Expt EFM PAGE 124 124 CHAPTER 8 DESIGN OF A MULTIBODY SYSTEM FOR REDUCED JOINT W EAR MAINTENANCE COSTS In Chapters 5 and 6, the procedure to an alyze multibody systems with joint wear was presented. The methodology in the analysis disc ussed in Chapter 5 was based on the contact force law and the finite element method (FEM), whereas the procedure pr esented in Chapter 6 was based on the elastic foundation model (EFM). These procedures were later validated against experiments in Chapter 7 and their performances ar e evaluated. In this Chapter, application of the FEM based procedure in the design of a mechanical system will be presented. Introduction Most mechanical/multibody systems consist of multiple joints that connect the system components. For example, the backhoe system shown in Figure 81 consists of three main revolu te joints which connect the bucket to the di pper, the dipper to the boom and the boom to the main body, respectively. For this particular system wear is exp ected to occur at all the joint. If one cycle is defined as scoopi ng and dropping off the dirt using the bucket then the wear depth at any of the joints is a function of 1) the num ber of cycles of operati on, 2) the joint reaction forces, 3) the material used to construct the joint com ponents (pin and bushing), 4) the size of the joint components (diameter and depth of the joint components) and 5) the operating conditions at the joint. In conventional design of such system s, it is unlikely that th e maximum allowed wear at all the joints will simultane ously occur at the same number of operation cycles. As a consequence, when the maximum allowable wear occurs on one of the jo ints the use of the system will have to be halted and the joint has to be repaired or replaced. The system can then be operated until the maximum allowable wear occurs on the other joint. The system operation is once again halted and the second joint is repaired or replaced. In the case of the backhoe system, repair of all the main joints w ould require that the system be out of operation in three occasions PAGE 125 125 and thus leading to high maintenance cost and lo ss of revenue (since th e equipment is out of service). Another alternative is that whenever one joint reaches the maximum allowable depth, the system operation is halted and all the joints are replace. For purposes of illustration it assumed that this is not the case. One design requirement could be to design the jo ints of the backhoe system such that the maximum allowable wear on all the joints occur at the same (or approximately the same) number of operation cycles. The result of this design is that the maintena nce cost and loss of revenue is reduced since the system will be out of operation only once as opposed to three times. This concept will be illustrated using a slidercra nk mechanism with wear occurring at two joints. Design Example: Design of a SliderCrank for Reduced Maintenance Cost The slidercrank mechanism used in Chapters 2, 3, 5 and 6 will be used to illustrate the design example. In this sliderc rank, however, wear is allowed to occur at two joint namely; the joint between the crank and the connectingrod which will be referred to as Joint 1 and the joint between the connectingrod and the slider referred to as Joint 2. A diagram of the slidercrank is shown in Figure 82 and its dimensions are shown in Table 81. Similar to the previous exam ples, it is assumed that the two joints consis t of a steel pin and PTFE bushing and thus wear on the pin is negligible when compared to the wear on the bushing. The material properties including the friction and wear coefficients for the joint components are listed in Table 82. Problem Definition The design problem can be stated as follows: given that the crank is operated at 30rpm, design the bushing (inner diameter and depth) for the two joints so that the maximum allowable wear at the two joints is attained simultaneousl y (at the same number of cycles of operation). This problem can be stated mathematically as follows: PAGE 126 126 2 1212111222 12 12Minimize : ,,, Such that : ,LB UB LB UBfddhhcycdhcycdh dddd hhhh (81) where cyc1 and cyc2 are the number of operation cycles required for wear to accumulate to the maximum allowable amount for Joint1 and Join t2, respectively, and the objective function f is the square of the difference between cyc1 and cyc2. In Eq. (81), the design variables are the dim ensions of the bushing; i.e., the diameter and the depth of the bushings ( d1 d2 h1 h2). The variables are constrained betw een their corresponding upper ( dUB hUB) and lower bounds (dLB hLB). The solution of the design problem will be deco mposed into two parts. In the first part, the analysis of the slidercrank with two joints we aring will be discussed. This will involve the integration of the dynamic analysis of the mechanis m and the wear analysis at the two joints. The second part will involve the solution of the optimization problem stated in Eq. (81). As will be seen la ter, the objective function is quite expe nsive to evaluate and thus a surrogate based optimization technique will be employed. Analysis of the SliderCrank with Multiple Joints Wearing The procedures presented for the analysis of multibody systems with joint wear consisted of two parts, namely; 1) dynamic analysis that accounts for joint clearance and 2) wear analysis to predict the amount of wear based on the preceding dynamic analysis. The dynamic analysis involves assembling the differential algebraic eq uations (DAE) of motion in which the system kinematic constraints and applie d loads should be specified. For this particular example the kinematic constraints can be expressed as shown in Eq. (82). In Eq.(82) the generalized coordinates are sim ilar to those described in Figure 36. It is worth mentioning that only the driving constraint, the cranks cons traint to the ground, and the rota tional constraint of the slider PAGE 127 127 appear in Eq. (82). The kinematic constraints that w ould norm ally have represented Joint1 and Joint2 (for perfect joints) are now replaced with force constraints described in Eq. (316). 111 111 3 3 1cos 0 sin 0 = 0 0 0 xl yl y t (82) The kinematic constraint and the force constrai nts can be assembled into the DAE which in turn can be solved to reveal the system dyna mics. Representative results from the dynamic analysis of this system are shown in Figure 83, Figure 84 and Figure 85. For purposes of illus tration and comparison, bushing diameters of 40mm and clearances of 0.0005mm was used (for both joints) to generate these results. In Figure 83 the reaction for c es for Joint1 and Joint2 have been plotted. The force response for Join t1 is identical to what was obtained for the dynamics based on the contact force model and the EFM based model (see Figure 66). This was expected since the clearances in the two joints are sm all enough so that the systems are also identical. In Figure 84, a plot of the incremental slidi ng angle is shown for both the joints. This is different f rom the incremental siding distan ce which would be a function of the bushing/pin radius. It should be pointed out that for a larger radius the incremental sliding distance would increase and so would the wear. The locus of the center of the contact region (for one cycle) is shown in Figure 85. For the same reasons, these re sults are once again identical to those observed in Figure 54 for the contact force based model and Figure 67 for the EFM based model. Solution of Optimization Problem Upon careful examination of Eq. (81), one will realize that the objective function is extrem ely expensive to evaluate. This is because determining a single value of cyc1 or cyc2 for PAGE 128 128 any combination of the joint diameters and depth ( d1 d2 h1 h2) requires a complete wear analysis involving several thousand cycles. Since numerous evaluation of the objective function is required during optimization, the current implicit form of the objective function will render the optimization problem unreasonably costly. Instead of conducting the optim ization using the current form of the objective function, also called the high fidelity model/function, a su rrogate based optimization approach [9196] can be employed. This approach entails constructing surrogate models, such as the response surface approximations [9799], support vector regressi on [100 101], and krigi ng [102 103], using data obtained from the highfidelity models. The su rrogates replace the highfidelity objective functions (and constraints) and thus offer fast approximations of the objective functions (and constraints) at other locations in the design space. The consequence is that the speed of optimization is substantially increased but at th e cost of accuracy. There are however procedures that can be used to improve the surrogate base d optimization results such as the trust region approach in which construction of the surrogate is focused in regions of possible optima based on previous optimization results [96 104]. The process of drawing data from the highfide lity models to construct the surrogate model is referred to as the design of experiments (DOE ). It should be mentioned that the choice of DOE will have a great influence on the quality of the su rrogate model. The reader is referred to the works by Simpson et al [105], Gi unta et al [106] and Goel [107] for a detailed discussion on the subject. Despite the various possible su rrogate models available, for this example, the objective function will be replaced by a response surface a pproximation. Furthermore the design points are selected by a space filling technique called Latin Hypercube Sampling (LHS) [108]. For PAGE 129 129 convenience the surrogate toolbox developed by Viana [109] was us ed to generate the DOE as well as the surrogate models. For the current slidercrank problem the maxi mum allowable wear depth for both joints was set at 2mm and the design range for the bus hing diameters (for the two joints) was set between 20mm and 70mm, whereas the range for the bushing depth was set between 10 mm and 50mm. These specifications are summarized in Table 83. The design space is defined as any com bination of the two variables with the two ra nges. A plot of the de sign points, generated using LHS is shown in Figure 86. In this plot the design space is boxlike w ith both variables (diam eter and depth) normalized between 0 and 1. It is emphasized that this DOE was used to construct the surr ogate models for cyc1 and cyc2 in Eq.(81). In Figure 87 and Figure 88, the response of cyc1 and cyc2 generated using the constructed surrogate models for cyc1 and cyc2 are shown. It can be inferred from both plots that as the bushing depth increases, cyc1 and cyc2 also increase. This is r easonable since a longer bushing depth would mean that the contact pressure re sponsible for wear is distributed over a longer length. As a result, the rate of cha nge of the wear depth would reduce and cyc1 and cyc2 would increase. It can also be in ferred from both figures that as the diameter increases, cyc1 and cyc2 decreases. This observation seems counterintuitive since a larger diameter would generally mean that the resulting contact pressure is distributed over a greater pe rimeter and thus dictating that rate of change of the wear depth would reduce causing cyc1 and cyc2 to increase. In reality however, an increase in the diameter will cause the incremental siding distance to increase which will in turn cause the rate of wear depth to increase and cyc1 and cyc2 to decrease. The later behavior is however not as strong as the former. It is therefore true to state that the bushing depth has a dominating effect on the value of cyc1 and cyc2. PAGE 130 130 With the surrogates for cyc1 and cyc2 available, the optimization problem in Eq. (81) can be solved. F or this problem a standard optimizer in the Matlab software was used to solve the problem. The results from the optimization is shown in Table 84. The combination of bushing dim ensions (diameter and depth) for both joints necessary to minimize the difference between cyc1 and cyc2 are shown in the table. The difference from the optimization is zero cycles which means that (according to the surrogate models) th e maximum allowable wear depth in both joints will be achieved simultaneously. It is however e xpected that the results are not dead accurate since they were obtained using surrogate m odels which provide approximates for the highfidelity model. Listed in Table 85, are the results generated u sing the highfidelity model for the combination of bushing dimension at the optimum solution. For the optimum configuration the highfidelity model predicted a difference of 4000 cycles (3.3%) between cyc1 and cyc2. The ideal condition would require that the difference be equal to zero Also the highfidelity model predicted a difference of 2165 cycles a nd 1835 cycles from the surrogate for cyc1 and cyc2 respectively. These differences are attributed to the approximate nature of the surrogate model. These results can however be improved by using the trust regi on approach [96,104] Summary and Concluding Remarks In this Chapter an illustra tion of the application of th e procedure to analyze multibody systems with joint wear was presented. The a pplication involved the design of a slidercrank mechanism so as to allow the maximum allowabl e wear depth at two joints to be attained simultaneously. This was restated as an optimi zation problem that would involve minimizing the difference in the number of cycles required to attain the maximum wear depth for both cycles. The solution was obtained by minimizing a surrogat e model that was constructed in place of the actual objective function. PAGE 131 131 The application demonstrates an example of how the procedures presented can be use to improve the design of planar multibody system with joint wear. Other examples include minimizing joint wear by changing the system mass and dimensions, reducing out of plane motion cause by the wear by redesigning the system predicting amount of wear at the joints of system for scheduling maintenance or issuing warra nties etc. The procedur es may therefore be a useful tool for the designer. PAGE 132 132 Table 81. Dimension and mass pa rameter for slidercrank mechanism Length (mm) Mass (g) Moment of inertia (kgm2) Crank 1.00 10.00 45.00 Connecting Rod 1.75 15.00 35.00 slider 30.00 Table 82. Material proper ties for the joint components Pin Bushing Youngs modulus 0.29 0.38 Poisson ratio 206.8 GPa 0.5 GPa Wear coefficient (pin & bushing[75] 5.05x104 mm3/Nm Friction coefficient (pin & bushing) [46] 0.13 Table 83. Parameter and design space specifications fo r the optimization Joint 1 Joint2 Maximum allowable wear depth 2 mm 2 mm Bushing diameter design range 20 mm d1 70 mm20 mm d2 70 mm Bushing depth design range 10 mm h1 50 mm10 mm h2 50 mm Table 84. Solution of optimization problem (Eq. (81)) Value Value of objective ( f ) 0 cycles Bushing diameter Joint 1 (1d ) 46.58 mm Bushing diameter Joint 2 (2d ) 48.60 mm Bushing depth Joint 1 (1h ) 37.25 mm Bushing depth Joint 2 (2h ) 14.00 mm Table 85. Comparison of results betw een surrogate and highfidelity model cyc1 cyc2 Difference Surrogate 119435119435 0 (0%) Highfidelity (actual wear simulation) 1216001176004000 (3.3%) Difference 2165 (1.8%)1835 (1.5%) PAGE 133 133 Figure 81. Backhoe system w ith three main revolute joints. Figure 82. Slidercrank mechanism with two wearing joints wearing. Crank Connecting rod Wearin g Joint 1 Groun d Ground Slide r Wearin g Joint 2 Joint 3 Joint 1 Joint 2 Di pp e r Boo m Bucket PAGE 134 134 0 1 2 3 4 5 6 0 1000 2000 3000 4000 5000 Crank Position (rad)Force (N)Joint Force for Joint 1 & Joint 2 Joint 1 Joint2 Figure 83. Joint reaction force for Joint1 and Joint 2 0 1 2 3 4 5 6 0 0.005 0.01 0.015 0.02 0.025 0.03 Crank Position (rad)Incremental Sliding Angle (rad)Incremental Sliding Angle Joint 1 Joint 2 Figure 84. Incremental sliding angle Joint 1 and Joint 2. PAGE 135 135 30 20 10 0 10 20 30 30 20 10 0 10 20 30 xcoordinate (mm)ycoordinate (mm)Locus of the Center of the Contact Region (Joint1) Pin Bushing 30 20 10 0 10 20 30 30 20 10 0 10 20 30 xcoordinate (mm)ycoor di na t e ( mm ) Locus of the Center of the Contact Region (Joint2) Pin Bushing Figure 85. Locus of the center of the contact region. a) Locus for Joint 1 b) Locus for Joint 2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 DiameterDepthLHS (Normalized Design Space) Figure 86. DOE used to construct surrogate for cyc1 and cyc2. PAGE 136 136 0 0.5 1 0 0.5 1 0 1 2 3 4 x 105 Diameter Number of Cycles to Reach Max. Allowable Wear (Joint 1) Depth Cycles Figure 87. Response genera ted using the surrogate for cyc1. 0 0.5 1 0 0.5 1 0 2 4 6 8 x 105 Diameter Number of Cycles to Reach Max. Allowable Wear (Joint 2) Depth Cycles Figure 88. Response genera ted using the surrogate for cyc2. PAGE 137 137 CHAPTER 9 SUMMARY AND FUTURE WORK In this Chap ter, a summary of the main ideas in this work will be presented. In addition, suggestions for future research will be discussed. Summary and Discussion Analysis of planar multibody systems with revo lute joint wear was the central theme of this work. Two procedures to analyzing such syst ems were developed. In the first procedure the analysis is based on a contact forces law, used to predict the joint reactio n forces, and the finite element method employed for the purpose of wear prediction. In the second procedure the elastic foundation model is used to determine the contact forces as well as for the wear prediction. The basic framework for theses procedures consists of three main parts. These are; 1) modeling a nonideal revolute joint that would replicate the behavior of a join t with clearance, 2) developing an efficient procedure to predict wear and 3) integrating the wear pr ocedure into the dynamic analysis of multibody systems. The three parts will briefly be summarized in the context of the two procedures. For the two procedures, a nonideal revolute jo int (planar) was modeled using a procedure that closely resembles a real nonideal joint. Unlike the ideal joint which uses kinematic constraints to restrict the motion of the jo int components, the nonid eal joint uses force constraints to guide the motion of the joint co mponents. The procedure assumes that the joint components can exhibit three possible configurati ons, namely: 1) freeflight, 2) impact and 3) following motion. In the case of the freeflight, no contact occurs between the components. This configuration is modeled by requiri ng that the contact/reaction for ce be zero. Thus no restriction is imposed on the motion of the joint components. For the impact and following motion, contact PAGE 138 138 between the joint components is established. Contact force is thus developed based on the amount of penetration experienced during contact. In the first procedure, the contact force is computed using a modified form of the Hertz contact law, whereas in the second procedure the contact force is determ ined using the elastic foundation model. In the latter cas e, the contact pressure between the joint components is first determined, and thereafter, the resultant contact force can be evaluated. In both procedures, however, the magnitude of the force is a functio n of penetration between the joint components. To enforce the joint constraint, the contact fo rce is applied on the components to prevent further penetration. Due to this contact force, the motion of the compone nts is restricted. The implementation of this joint was illustrated with the aid of a slidercrank mechanism with a nonideal revolute joint between the crank and the connecting rod. Using the mechanism, it was shown that the system dynamics is altered as the joint clearance size vari es. The two procedures showed reasonable agreement in the force magn itude for a small clearance but appreciable disagreement in the magnitude for relatively large clearances. However, the trend in the force distribution for the force dist ribution remained the same. A methodology to predict wear that is built upon a widely used iterative prediction procedure was discussed. In the iterative proce dure, the wear occurring at the contact interface between two bodies that are in relative motion is determined based on Archards wear model. Incremental wear is determined at each iterati on and accumulated up to the desired number of cycles. In every iteration (or several iterations) the geometry is updated to reflect the worn material. In the first procedure (for analysis of multibody systems with joint wear) the finite element method is used to determine the contact pressure necessary fo r wear calculation, whereas in the second procedure th e contact pressure is readily available from the analysis used PAGE 139 139 to determine the contact force. In order to redu ce computational costs, an adaptive extrapolation scheme is used. The two main contributions in the wear prediction pr ocedures include: 1) development of an adaptive extrapolation scheme in which a systematic approach in selecting the extrapolation was proposed and 2) the developm ent of a updating procedure that maintains smooth boundary and ensures regularity of the fi nite element mesh at the boundary (in the case of the FEM based procedure). Finally the analysis of multibody systems with joint wear is completed by integrating the wear prediction procedures into the dynamic anal ysis. In the integration process, the dynamic analysis provides information about the system dynamic for multibody systems with joint clearance. This information includes, contact force (and contact pressure for the second procedure), location of contact a nd the incremental siding distance. On the other hand, the wear prediction determines the amount of wear at the joints based on th e result of the system dynamics and the contact geometry is updated after every cycle to reflect the wo rn out material. The changes in the geometry due to the wear are accounted for in the dynamic analysis by updating the clearance. In this case the clearance is a vector that is dependent on the bushing angular coordinate. Thus any change in the system dynami cs (contact force and c ontact locations) due to the wear is captured in the integration process. The use of the integrated model for the two procedures was demonstrated using a slidercrank mechanism in which the joint between th e crank and the connecting rod was allowed to wear. It was found that the two procedures predicted different wear profiles and different maximum wear depth but the same location of the maximum wear and the same wear volume. The difference in the wear profiles is attributed to the simplification in the EFM where contact pressure is determined without considering how displacements at one co ntact region affects the PAGE 140 140 displacement at other regions (ind ependent springs). On the other hand the equality of the wear volume for the two models is a manifestation of the equality in the force predicted by the two models. Experiments were conducted to validate the tw o procedures and to assess the performance between the two procedures. An experimental slidercrank mechanism was used to facilitate the validation. The validation revealed that both procedures predicted reasonably accurate wear volumes and location of maximum wear. With re gard to the wear profile and maximum wear depth, the FEM based procedure gave more accurate predictions than the EFM based procedure. On the other hand the EFM based procedure had a much faster computation time than the FEM based procedure. The conclusion that can be drawn from the validation is that the FEM procedure will be a better procedure in the anal ysis of multibody systems with joint wear when the computational costs in not of concern. On th e other hand if the cost is of concern and only qualitative information about the system is n eeded then the EFM based procedure will be the suitable choice. Other scenarios will require a co mpromise on either accuracy or computational costs. Future Work This work has primarily focused on the analys is of systems for which the joints are nonlubricated. When considering real systems, this condition can be considered to be among the extreme cases since in most real systems some t ype of lubrication is used. Thus, this works will serve as a starting or as a refere nce point for other realistic cases. Naturally, the lubricated joint should be the next case considered. In this case the lubricant is expected to affect the system dynamics. The dynamic analysis of multibody system s with joint lubricant has already been explored [14 15 18, 111,112]. It is also expected that the lubricant wi ll affect the wear mechanism and thus the wear rate. PAGE 141 141 Another area interest that will follow on from the lubricated joint is the lubricated joint with trapped impurities. This case is probably the most realistic case of all three mentioned. In most joint encountered, the join t is lubricated. However, while the system is in operation, impurities such as sand may become trapped betw een the joint components together with the lubricants. The presence of the im purities in the joint will most certainly have an effect on the wear and possible an effect on th e system dynamics (could serve to increase the friction or act as an additional lubricant). The final suggestion for future research is the analysis of spatial multibody systems with joint wear. This can be coupled with joint lubri cants with impurities to ma ke it the most general and realistic case. PAGE 142 142 LIST OF REFERENCES 1. Dubowsky, S., Freudenstein, F ., 1971, Dynamic Analysis of Mechanical Systems with Clearances Part1: Formulation of Dynamic mo del, Journal of Engineering for Industry, pp. 305. 2. Dubowsky, S., Freudenstein, F ., 1971, Dynamic Analysis of Mechanical Systems With Clearances, Part 2: Dynamic Response, Journal of Engi neering for Industry, pp. 310 316. 3. Dubowsky, S., Gardner, T.N., 1977, Desi gn and Analysis of Multilink Flexible Mechanisms with Multiple Clearance Conn ection, Journal of Engineering, pp. 88. 4. Earles, S.W.E., Wu, C.L.S., 1973, Motion An alysis of a Rigid Link Mechanism with Clearance at a Bearing Using Lagrangian Mechanics and Digital Computation Mechanisms, pp. 83. 5. Wu, C.L.S., and Earles, S.W.E ., 1977, A Determination of ContactLoss at a Bearing of a Linkage Mechanism, Journal of Engineering for Industry, Series B, 99(2), pp. 375 380. 6. Furuhashi, T., Morita, N., Matsuura, M., 1978, Research on Dynamics of FourBar Linkage with Clearances at Turning Pairs 1st Report: General Theo ry Using Continuous Contact Model, JSME, 21, pp. 518. 7. Morita, N., Furuhashi, T., Matsuura, M., 1978, Research on Dynamics of FourBar Linkage with Clearances at Turning Pairs 2nd Report: Analysis of CrankLever Mechanism with Clearance at Joint of Crank Coupler Usi ng Continuous Contact Model, JSME, 21, pp. 1284. 8. Morita, N., Furuhashi, T., Matsuura, M., 1978, Research on Dynamics of FourBar Linkage with Clearances at Turning Pairs 3rd Report: Analysis of CrankLever Mechanism with Clearance at Joint of Coupler and Lever Using Continuous Contact Model, JSME, 21, pp. 1292. 9. Furuhashi, T., Morita, N., Matsuura, M., 1978, Research on Dynamics of FourBar Linkage with Clearances at Turning Pairs 4th Report: Forces Acting at Joints of CrankLever Mechanism, JSME, 21, pp. 1299. 10. Farahanchi, F., Shaw, S. W., 1994, Chaotic and Periodic Dynamics of a SliderCrank Mechanism with Slider Clearance, Journal of Sound and Vibration, 177(3), pp. 307 324. 11. Rhee, J., Akay A., 1996, Dynamic Response of a Revolute Joint with Clearance, Mechanisms Machine Theory, 31(1), pp. 121. PAGE 143 143 12. Khulief, Y.A., and Shabana, A.A., 1986, Dyna mic Analysis of Constrained System of Rigid and Flexible Bodies with Intermittent Motion, Journal of Mechanisms, Transmissions, and Automations in Design, 108, pp.38. 13. Ravn, P., 1998, A Continuous Analysis Me thod for Planar Multibody Systems with Joint Clearance, Multibody Systems Dyna mics, Kluwer Academic Publishers, 2, pp. 1 24. 14. Flores, P., Ambrsio, J., Claro, P.J., 2004, Dynamic Analysis for Planar Multibody Mechanical Systems with Lubricated Joints, Multibody System Dynamics, 12, pp. 47 74. 15. Ravn, P., Shivaswamy, S., Alshaer, B. J., La nkarani, H. M., 2000, Joint Clearances With Lubricated Long Bearings in Multibody Mechanical System s, Journal of Mechanical Design, 122, pp. 484. 16. Flores, P., Ambrsio, J., 2004, Revolute Join ts with Clearance in Multibody Systems, Computers and Structures, 82, pp. 1359. 17. Flores P., 2004, Dynamic Analysis of Mechanical Systems with Imperfect Kinematic Joints, Ph.D. Thesis, Minho Univ ersity (Portugal), Guimares. 18. Flores, P., Ambrsio, J., Claro, J. C. P., Lankarani, H. M.Koshy, C. S., 2006, A Study on Dynamics of Mechanical Systems Includi ng Joints with Clearance and Lubrication, Mechanism and Machine Theory, 41, pp. 247. 19. Flores, P., 2009, Modeling and simulation of wear in revolute clearance joints in multibody systems, Mechanism and machine theory, 44(6), pp. 1211. 20. Nikravesh, P.E., 1988, ComputerAided Analys is of Mechanical System, PrenticeHall, Englewood Cliffs, NJ. 21. Haug, E. R., 1989, ComputerAided Kinematic s and Dynamics of Mechanical Systems, Allyn and Bacon, Needham Heights, USMA. 22. Shabana, A.A., 1989, Dynamics of Mu ltibody Systems, Wiley, New York. 23. Wittenburg, J., 2007, Dynamics of Multibody Systems, Springer. 24. Schiehlen, W., 1997, Multibody System Dynamics: Roots and Perspectives, Multibody System Dynamics, 1, pp. 149. 25. Atkinson, K.E., 1978, An Introduction to Numerical Analysis, Wiley, New York. 26. Baumgarte, J., 1972, Stabilization of Constr aints and Integrals of Motion in Dynamical Systems, Computer Methods in App lied Mechanics and Engineering, 1, pp. 1. PAGE 144 144 27. Lin, S.T., Huang J.N., 2002, Stabilizati on of Baumgartes Method Using the RungeKutta Approach, Journal of Mechanical Design, 124, pp. 633. 28. Lin, S.T., Huang J.N., 2002, Numerical In tegration of Multibody Mechanical Systems Using Baumgartes Constraint Stabilization Method, Journal of the Chinese Institute of Engineers, 25, pp. 242. 29. Wehage, R.A., Haug, E.J., 1982, Generalized Coordinate Partitioning for Dimension Reduction in Analysis of Constrained Dynami c Systems, Journal of Mechanical Design, 104(1), pp. 247. 30. Park, T., 1986 A Hybrid Constraint Stabiliz ationGeneralized Coor dinate Partitioning Method for Machine Dynamics, Journal of Mechanisms, Transmission, and Automation in Design, 108(2), pp. 211. 31. Bayo, E., Garcia de Jalon, J., and Serna, A., 1988, A Modified Lagrangian Formulation for the Dynamic Analysis of Constrained M echanical Systems, Computer Methods in Applied Mechanics and Engineering, 71, pp. 183. 32. Dubowsky, S., 1974, On Predicting the Dyna mic Effects of Clearances in Planar Mechanisms, ASME, Journal of E ngineering for Industry, pp. 317. 33. Ting, K.L., Zhu, J., Watkins, D., 2000The effects of joint clearance on position and orientation deviation of li nkages and manipulators, Mechanism and Machine Theory, 35(3), pp 391. 34. Kakizaki, T., Deck, J.F., Dubowsky, S., 1993, Modeling the Spatial Dynamics of Robotic Manipulators with Fixable Links a nd Joint Clearances, ASME, Journal of Mechanical Design, 115, pp. 839. 35. Farahanchi, F., Shaw, S. W., 1994, Chaotic and Periodic Dynamics of a SliderCrank Mechanism with Slider Clearance, Journal of Sound and Vibration, 177(3), pp. 307 324. 36. Schwab, A. L., Meijaard, J.P., Meijers, P., 2002, A Comparison of Revolute Joint Clearance Models in the Dynamic Analysis of Rigid and Elastic Mechanical Systems, Mechanisms and Machine Theory, 37, pp. 895. 37. Bauchau, O. A., Rodriguez, J., 2002, Modeli ng of Joints with Clearance in Flexible Multibody Systems, Journal of Solids and Structures, 39, pp 41. 38. Tsai, M.J., Lai, T.H., 2004 Kinematic Sensitivity Analysis of Linkage with Joint Clearance Based on Transmission Quality, Mechanism and Machine Theory, 39, pp. 1189, 39. Garcia Orden, J. C., 2005, Analysis of Joint Clearances in Multibody Systems, Multibody System Dynamics, 13, pp. 401. PAGE 145 145 40. Mukras, S., Mauntler, N., Kim, N.H., Schmitz, T.L., Sawyer, W.G., 2008, Dynamic Modeling of a SliderCrank Mechanism w ith Joint Wear, Paper No. DETC200849244, ASME 32nd Annual Mechanism and Robotic s Conference, New York, New York. 41. Alessandro, T., Edzeario, P., Marco, S., 2004, Experimental Investigation of Clearance Effect in A Revolute Joint, AIMETA Intern ational Tribology Conf erence, Rome, Italy. 42. Chunmei, J., Yang, Q., Ling, F., Ling, Z., 2002, The NonLinear Dynamic Behavior of an Elastic Linkage Mechanism with Clear ances, Journal of Sound and Vibration, 249(2), pp. 213. 43. Lankarani, H.M., Nikravesh, P.E., 1994, Cont inuous Contact Force Models for Impact Analysis in Multibody Systems, Nonlinear Dynamics, 5, pp. 193. 44. Khulief, Y.A., Shabana, A.A., 1987, A continuous force model for the impact analysis of flexible multibody system s, Mechanism and Machine Theory, 22(3), pp. 213. 45. Hunt, K.H., Crossley, F.R., 1975, Coefficient of restitution interpreted as damping in vibroimpact, Journal of Applied Mechanics, 7, pp. 440. 46. Schmitz, T.L., Action, J.E., Zigert, J.C ., Sawyer, W.G., 2005, The Difficulty of Measuring Low Friction: Uncertainty Analysis for Friction Coefficient Measurment, Journal of Tribology, 127, pp. 673. 47. Stolarski, T.A., 1990, Tribology in Machine Design, Heinemann Newnes, Jordan Hill, Oxford. 48. Sawyer, W.G., 2001, Wear Predictions for a SimpleCam Including the Coupled Evolution of Wear and Load, Lubrication Engineering, pp. 31. 49. Blanchet, T.A., 1997, The Inte raction of Wear and Dynamics of a Simple Mechanism, Journal of Tribology, 119, pp. 597. 50. Kim, N.H., Won, D., Burris, D., Holtkamp, B., Gessel, G.R ., Swanson, P., Sawyer, W.G., 2005, Finite element analysis and experiments of metal/metal wear in oscillatory contacts, Wear, 258, pp 1787. 51. Oqvist, M., 2001, Numerical simulations of mild wear using updated geometry with different step size approaches, Wear, 249, pp. 6. 52. Podra, P., and Andersson, S., 1997, Wear simulation with the Winkler surface model, Wear, 207, pp. 79. 53. Podra, P., and Andersson, S., 1999, Finite el ement analysis wear simulation of a conical spinning contact considering surface topography, Wear, 224, pp. 13. 54. Podra, P., and Andersson, S., 1999, Simul ating sliding wear with finite element method, Tribology International, 32, pp. 71l. PAGE 146 146 55. Mukras, S., Kim, N.H., Sawyer, W.G., Jackson, D.B., Bergquist, L.W., 2009, Numerical integration schemes and parallel computati on for wear prediction using finite element method, Wear, 266 pp. 822. 56. Flodin, A., and Andersson, S., 1997, Simula tion of Mild Wear in Spur Gears, Wear, 207, pp. 16. 57. Flodin, A., and Andersson, S., 2001, A simp lified model for wear prediction in helical gears, Wear, 249, pp 285. 58. Brauer, J., and Andersson, S., 2003, Sim ulation of Wear in Gears with Flank Interferencea Mixed FE and Analytical Approach, Wear, 254, pp 1216. 59. Hugnell, A., and Andersson, S., 1994, Simulating follower wear in a camfollower contact, Wear, 179, pp. 101. 60. Hugnell, A.B.J., Bjrklund, S., and Andersson, S., 1996, Simulation of the Mild Wear in a CamFollower Contact with Follower Rotation, Wear, 199, pp. 202. 61. Nayak, N., Lakshminarayanan, P.A., Babu, M.K.G., and Dani, A.D., 2006, Predictions of cam follower wear in diesel engines, Wear, 260, pp. 181. 62. Fregly, B.J., Sawyer, W.G., Harman, M. K., and S. Banks, A., 2005, Computational wear prediction of a total knee replacemen t from in vivo kinematics, Journal of Biomechanics, 38, pp. 305. 63. Maxian,T.A., Brown, T.D., Pedersen, D.R., and Callaghan, J.J., 1996,A SlidingDistance Coupled Finite Element Formula tion for polyethylene wear in total hip arthroplasty, Journal of Biomechanics, 29, pp. 687. 64. Bevill, S.L., Bevill, G.R., Penmetsa, J.R., Petrella, A.J., and Rullkoetter, P.J., 2005, Finite Element Simulation of Early Creep and Wear in Total Hip Arthroplasty, Journal of Biomechanics, 38, pp. 2365. 65. McColl, I.R., Ding, J., and Leen, S.B., 2004, Finite element simulation and experimental validation of fretting wear, Wear, 256, pp. 1114. 66. Dickrell, D.J., Dooner, III,D.B., and Sawyer, W.G., 2003, The evolution of geometry for a wearing circular cam: analytical and computer simulation with comparison to experiment, ASME Journal of Tribology 125, pp. 187. 67. Yan, W., O'Dowd, N.P., and BussoE.P., 2002, Numerical study of slide wear caused by a loaded pin on a rotating disc, Journa l of Mechanics and Physics of Solids, 50, pp. 449. 68. Archard, J.F., 1953, Contact and rubbing of flat surfaces, J Appl. Phys. Vol. 24, pp. 981. PAGE 147 147 69. Holm, R., 1946 Electric Contacts Uppsala: Almqvist & Wiksells Boktryckeri. 70. Cantizano, A., Carnicero, A., Zavarise, G., 2002, Numerical simulation of wearmechanism maps, Computat ional Materials Science, 25 pp.54. 71. Hegadekatte, V., Huber, N., Kraft, O., 2005, Finite element based simulation of dry sliding wear. IOP Publishing, 13, pp. 57. 72. Gonzalez, C., Martin, A., Llorca, J., Garrido, M.A., Gomez, M.T., Rico, A., Rodriguez, J., 2005, Numerical Analysis of PinonDi sk Test on AlLi/SiC composites, Wear, 259, pp. 609. 73. Sfantos, G.K., Aliabadi, M.H., 2006, Wear simulation using an incremental sliding boundary element method, Wear 260, pp. 1119. 74. Telliskivi, T., 2004, Simulation of wear in a rollingsliding contact by a semiWinkler model and the Archards wear law, Wear 256, pp. 817. 75. Schmitz T.L., Action, J.E., Burris, D.L., Zi egert, J.C., Sawyer, W.G., 2004, WearRate Uncertainty Analysis, AS ME Journal of Tribology, 126(4), pp. 802. 76. Yang, L.J, 2005, A test methodology for dete rmination of wear coefficient, Wear, 259 pp. 1453. 77. Sawyer, W. G., 2004, Surface shape and cont act pressure evolution in two component surfaces: application to copper chemicalm echanicalpolishing, Tribology Letters, 17(2) pp. 139. 78. Choi, K.K., Chang, K.H., 1994, A study of design velocity field computation for shape optimal design, Finite Elements in Analysis and Design, 15(4), pp. 317. 79. Bei, Y., and Fregly, B.J., 2004, Multibody dynamic simulation of knee contact mechanics, Medical Engineering and Physics 26, pp. 777. 80. Johnson, K.L., 1985, Contact mechanics Cambridge University Press, Cambridge. 81. Fregly, B.J., Bei, Y., and Sylvester, M.E., 2003, Experimental evaluation of an elastic foundation to predict contact pressure in knee replacement, Journal of Biomechanics 36(11) pp. 1659. 82. Blankevoort, L., Kuiper, J.H., Huiskes, R., Gr ootenboer, H.J., 1991, Articular contact in a threedimensional model of the knee, Journal of Biomechanics, 24, pp. 1019. 83. Pandy, M.G., Sasaki, K., Kim, S., 1997, A threedimensional musculoskeletal model of the human knee joint. Part 1: theoreti cal construction, Computer Methods in Biomechanics and Biomedical Engineering, 1, pp. 87. PAGE 148 148 84. Pandy, M.G., Sasaki K., 1998, A threedimensional musculoskeletal model of the human knee joint. Part 2: analysis of ligament function, Computer Methods in Biomechanics and Biomedical Engineering, 1, pp. 265. 85. An, K.N., Himenso, S., Tsumura, H., Kawai, T., Chao, E.Y.S., 1990, Pressure distribution on articular surfaces: application to joint stability analysis, Journal of Biomechanics, 23, pp. 1013. 86. Li, G., Sakamoto, M., Chao, E.Y.S., 1997, A comparison of different methods in predicting static pressure distribution in articulating joints, Journal of Biomechanics, 30, pp. 635. 87. Fregly, B.J., Sawyer, W.G., Harman, M.K., Banks, A.S., 2005, Computational wear prediction of a total knee replacement from in vivo kinematics, Journal of biomechanics, 38(2), pp. 305. 88. Halloran, J.P., Easley, S.K., Petrella, A. J., Rullkoetter, P.J., 2005 Comparison of deformable and Elastic foundation finite element simulations for predicting knee replacement mechanics Journal of Biomechanical Engineering 127 813 818 89. Mauntler, N., Kim, N.H., Sawyer, W.G., Schmitz, T.L., 2007, An instrumented crankslide r mechanism for validation of a comb ined finite element and wear model, Proceedings of 22nd Annual Meeting of Amer ican Society of Precision Engineering, October 14, Dallas, Texas. 90. Mauntler, N., Mukras, S., Kim, N.H., Schmitz, T. L., Sawyer, W. G., 2009, An Instrumented CrankSlider Mechanism for Wear Testing, STLE 64th Annual Meeting, Lake Buena Vista, Florida 91. Queipo, N.V., Haftka, R.T., Shyy,W., Goel,T ., Vaidyanathan, R., Tucker, P.K., 2005, Surrogatebased analysis and optimization, Progress in Aerospace Science, 41, pp 1 28. 92. Mack, Y., Goel, T., Shyy, W., Haftka, R.T., Queipo, N.V., 2005, Multiple surrogates for the shape optimization of bluff bodyfacilita ted mixing, Proceedings of the 43rd AIAA aerospace sciences meeting and exhibit, Reno, NV. Paper no., AIAA. 93. Booker, A.J., Dennis Jr, J.E., Frank, P.D., Serafini, D.B., Torczon V. and Trosset, M.W., 1999, A rigorous framework for optimizati on of expensive functions by surrogates, Structural Optimization, 17(1), pp. 1. 94. Audet, C., Dennis, Jr.,J.E., Moore, D.W., Booker, A.J, Frank ,P.D., 2000, A surrogatemodelbased m ethod for constrained op timization, 8th AIAA/USAF/NASA/ISSMO Symposium on Multidisiplanary Analysis a nd Optimization, Long Beach, CA, Paper no, AIAA20004891. PAGE 149 149 95. Jansson, T., Nilsson, L., Redhe, M., 2003, U sing surrogate models and response surfaces in structural optimization with applicati on to crashworthiness design and sheet metal forming, Structural Multidisciplinary Optimization, 25, pp. 129. 96. Mack, Y., Goel, T., Shyy, W., Haftka, R., 2007, Surrogate ModelBased Optimization Framework: A Case Study in Aerospace Design Studies in Computational Intelligence 51, pp. 323. 97. Roux, W.J., Stander, N., Haftka, R.T., 1996, Response surface approximations for structural optimization, AIAA/NASA/ISSMO, Symposium on Multidisciplinary Analysis and Optimization, Reston,VA. 98. Venter, G., Haftka, R.T., St arnes Jr, J.H., 1998Construc tion of Response Surface Approximations for Design Optimization, AIAA 36(12), pp. 2242. 99. Papila, M., Haftka R.T., 2000, Response surf ace approximations: noi se, error, repair and modeling errors, AIAA 38(12), pp. 2336 100. Smola, A.J., Schlkopf, B., 2004, A tutorial on support vector regre ssion, Statistics and computing, 14, pp. 199. 101. Girosi, F., 1998, An equiva lence between sparse approximation and support vector machines, Neural Computing, 10, pp.1455. 102. Martin, J.D., Simpson, T.W., 2004, On the use of kriging models to approximate deterministic computer models, ASME IDET Conferences and Computers and Information in Engineering Conference Salt Lake City, Utah. 103. Simpson, T.W., Mauery, T.M., Korte, J.J., Mistree, F., 1998, Comparison of response surface and kriging models for multidisciplin ary design optimizati on, AIAA, pp. 1. 104. Alexandrov, N.M., 1996, A trust region fram ework for managing approximation models in engineering optimization, AIAA/NAS A/USAF multidisciplinary analysis & optimization symposium, Bellevue, WA,. Paper no. 96. 105. Simpson, T.W., Lin, D.K.J., Chen, W., 2001, Sampling strategies for computer experiments: design and analysis International Journal of Reliability and Applications; 2(3):209. 106. Giunta, A.A., Wojtkiewicz, S.F., Eldred, M. S., 2003, Overview of modern design of experiments methods for computational simulations, Fortyfirst AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, 6 January; AIAA20030649. 107. Goel, T., Haftka, R.T., Shyy, W., Watson, L.T., 2008, Pitfalls of using a single criterion for selecting experimental designs, Int. J. Numer. Meth. Engng; 75:127. PAGE 150 150 108. Mckay, M.D., Beckman, R.J., Conover, W.J., 1979 A comparison of three methods for selecting values of input vari ables in the analysis of out put from a computer code, Technometrics, 21,pp. 239. 109. Viana, F.A.C., 2009, "Surrogates Toolbox Users Guide," Users manual, http://fchegury.googlepages.com. 110. Choi, K.K., Chang, K.H., 1994, A study of design velocity field computation for shape optimal design, Finite Elements in Analysis and Design, 15(4), pp. 317. 111. Alshaer, B.J., Nagarajan, H., Beheshti, H.K, Lankarani, H.M., Shivaswamy, S., 2005, Dynamics of a multibody mechanic al system with lubricated long journal bearings, J. Mech. Des., 127(3), 493. 112. Stefanelli, R., Valentini, P.P., Vita, L., 2005, Modeling of hydrodynamic journal bearing in spatial multibody system, Proceedings of IDETC/CIE, long beach, Ca. PAGE 151 151 BIOGRAPHICAL SKETCH Saad Mukras was born in Nairobi, K enya. He was raised in Nairobi and partially in Gaborone, Botswana, where he completed his secondary education. He then joined University of Botswana and then transferred to Embry Riddl e Aeronautical University in Daytona Beach, Florida. There, he studied airc raft engineering technology and rece ived his bachelors degree in 2003. In 2004 he joined the University of Florid a to pursue a masters and later a doctorate degree in mechanical engineer ing. He worked under the supervision of Dr. NamHo Kim, completing several research projects, earning hi s masters degree in 2006 and doctorate in 2009. 