• Author / Uploaded
• Bruno Ramos

Citation preview

Abstract LOWNDES, ERIK MCKENZIE.

Development of an Intermediate DOF Vehicle

Dynamics Model for Optimal Design Studies. (Under the direction of Dr. Joseph W. David)

The demands imposed by the optimal design process form a unique set of criteria for the development of a computational model for vehicle simulation. Due to the large number of simulations that must be performed to obtain an optimized design the model must be computationally efficient. A competing criterion is that the computational model must realistically model the vehicle. Current trends in vehicle simulation codes have tackled the problem of realism by constructing elaborate full vehicle models containing dozens if not hundreds of distinct bodies. Each body in a model of this type is associated with six degrees of freedom. Numerous constraint equations are applied to the bodies to represent the physical connections. While the formulation of the equations is not particularly difficult, and in fact has been automated in several software packages, the resulting model requires a considerable amount of computational time to run. This makes the model unsuitable for the application of computational optimal design techniques. Past research in the field of vehicle dynamics has produced numerous computational models which are small enough and fast enough to satisfy the speed

demands of the optimal design process. These models typically use less than a dozen degrees of freedom to model the vehicle. They do a good job of predicting the general motion of the vehicle and they are useful as design tools but they lack the accuracy required for optimal design. A model that bridges the gap between these two existing classes of models and is suitable for performing optimal design was developed. The model possesses twenty-eight degrees of freedom and consists of eight bodies which represent the sprung mass, the rear suspension, the left front spindle, the right front spindle, and the four wheels. A driver control algorithm was developed which is capable of driving the car near its handling limits. The NCSU Legends race car was modeled and an attempt was made to optimize the vehicle setup for the Kenley, NC race track.

Development of an Intermediate DOF Vehicle Dynamics Model for Optimal Design Studies

by

Erik M. Lowndes

A thesis submitted to the Graduate Faculty of North Carolina State University in partial fulfillment of the requirements for the Degree of

Doctor of Philosophy Department of Mechanical and Aerospace Engineering Raleigh 1998

Approved By:

for my wife, Celeste, and my son, Mason

ii

Biography The author was born in Boulder, Colorado on March 21, 1968. He graduated first from Oak Ridge High School, Oak Ridge, TN in 1986. He attended the University of Illinois at Urbana-Champaign in 1986 as a Chancellor’s Scholar in the Campus Honor’s Program. In January of 1990 he received his Bachelor of Science in Engineering Physics. He continued his studies at the University of Illinois, receiving a Master of Science in Physics in May of 1991. The author was married on June 23, 1991. In the Summer of 1992 the author entered the Ph.D. program in the Department of Mechanical Engineering at North Carolina State University, Raleigh, North Carolina. In March of 1994 he was elected to the honor society of Phi Kappa Phi. His first son, Mason, was born on April 23, 1997.

iii

Acknowledgments The author wishes to express his sincere appreciation to the chairman of his advisory committee, Dr. Joseph W. David, for his patience, encouragement, advise and support throughout this research. The author is also grateful to Dr. C. Tim Kelley, Dr. Larry M. Silverberg and Dr. John S. Strenkowski for serving on the advisory committee. Many thanks are also due to Dave Lewandowski and Mark Strohmeyer for their support, encouragement, and ideas. Additional thanks are due to Mark Etheridge and Jarno Kilian, who, in addition to their encouragement and valued input, loaned their personal computers to the author for use as part of the PVM parallel processing network used to generate the results in this thesis. Finally, the author would like to express his gratitude to his family, and especially to his wife Celeste and his son Mason, for their patient understanding, sacrifice and support (financial and otherwise) during the course of this study.

iv

Table of Contents List of Tables

viii

List of Figures

x

1

Introduction 1 1.1 Motivation for the Study ................................................................................1 1.2 Historical Background....................................................................................3 Vehicle Modeling ...........................................................................................6 Driver Modeling........................................................................................... 29 Model Parameter Measurement .................................................................... 39 Direct Measurement ............................................................................ 40 Parameter Identification....................................................................... 41 1.3 Derivation Methodology and Overview of the Thesis.................................... 43

2

Equations of Motion - Sprung Mass 47 2.1 Introduction ................................................................................................. 47 2.2 Sprung Mass Kinetic and Potential Energy Terms......................................... 48 2.3 Euler Parameter Constraints ......................................................................... 54

3

Equations of Motion - Front Suspension and Wheels 56 3.1 Introduction ................................................................................................. 56 3.2 Front Spindle Kinetic and Potential Energy Terms ........................................ 59 3.3 Front Wheel and Tire Rotational Energy Terms ............................................ 60 3.4 Generalized Forces for Springs and Dampers................................................ 62 3.5 Constraint Forces for the Control Arms ........................................................ 69 3.6 Summary of Results...................................................................................... 73

v

4

Equations of Motion - Three Link Rear Suspension 77 4.1 Introduction ................................................................................................. 77 4.2 Unsprung Mass Kinetic and Potential Energy Terms..................................... 77 4.3 Rear Wheel Rotational Energy Terms ........................................................... 79 4.4 Rear Springs and Dampers ........................................................................... 82 4.5 Panhard Rod and Trailing Link Constraints................................................... 84 4.6 Summary of Results...................................................................................... 86

5

Equations of Motion - Steering System 89 5.1 Introduction ................................................................................................. 89 5.2 Rack and Pinion ........................................................................................... 89 5.3 Four Bar Linkage ......................................................................................... 90 5.4 Tie Rod Constraints...................................................................................... 96

6

Road Model 99 6.1 Introduction ................................................................................................. 99 6.2 Road Surface Coordinate System ............................................................... 100 6.3 Location of the Tire to Road Contact Point ................................................ 102 6.4 Velocity of the Tire to Road Contact Point................................................. 108 6.5 Vehicle Position and Heading Angle ........................................................... 110 6.6 Road Segment Models................................................................................ 112 Linear Polynomial Road Segment ............................................................... 113 Quadratic Polynomial Road Segment.......................................................... 114 Cubic Polynomial Road Segment ................................................................ 115

7

Equations of Motion - Tire Model 117 7.1 Introduction ............................................................................................... 117 7.2 Coordinate Systems.................................................................................... 117 7.3 The Magic Formula Tire Model.................................................................. 119 Support Forces........................................................................................... 120 Tractive Forces .......................................................................................... 121 Lateral Slip and Longitudinal Slip ...................................................... 122 Magic Formula .................................................................................. 125 7.4 Generalized Force and Moments................................................................. 127 Front Tires ................................................................................................. 127 Rear Tires .................................................................................................. 129

vi

8

Equations of Motion - Driver Model 131 8.1 Introduction ............................................................................................... 131 8.2 Steering Control......................................................................................... 132 Driver Path Definition ................................................................................ 133 Steering Profile Optimization and Cost Function Computation.................... 134 8.3 Speed Control ............................................................................................ 138 Driver Dynamics Block .............................................................................. 140 Vehicle Dynamics Block............................................................................. 143 Preview Compensation Block ..................................................................... 144

9

Results, Conclusions and Recommendations 150 9.1 Introduction ............................................................................................... 150 9.2 Measurement Process and Model Data ....................................................... 150 Vehicle Data............................................................................................... 151 Tire Data.................................................................................................... 156 Track Data ................................................................................................. 160 Driver Model Data ..................................................................................... 161 9.3 Model Chassis Setup .................................................................................. 162 9.4 Simulation Results...................................................................................... 165 Optimizer and Cost Function Computation Setup ....................................... 165 Optimal Steering Profile Configuration ....................................................... 167 Optimal Velocity Profile Setup ................................................................... 168 Speed Control Algorithm Performance ....................................................... 170 Steering Control Performance..................................................................... 174 9.5 Vehicle Optimization Results...................................................................... 177 9.6 Recommendations for Future Research....................................................... 179

Bibliography

183

Appendix A - Useful Derivatives 190 Angular Velocity Derivatives .............................................................................. 190 Transformation Matrix Derivatives...................................................................... 192 Appendix B - Wheel Inertia Estimate 196 Thin Cylindrical Disk .......................................................................................... 196 Thin Walled Cylindrical Shell .............................................................................. 197 Rotating Assembly Model ................................................................................... 198 Appendix C - Tire Data 199 BFGoodrich Letter.............................................................................................. 199 Tire Data Plots.................................................................................................... 200

vii

List of Tables 1.1

Vehicle Model Degrees of Freedom ........................................................................ 2

1.2

Computational Degrees of Freedom ........................................................................ 3

1.3

Identification of the Sub-Terms in the Equations of Motion................................... 45

1.4

Organization of the Vehicle Model Derivations ..................................................... 46

3.1

Front Suspension Kinetic and Potential Energy Terms........................................... 73

3.2

Wheel and Tire Rotational Energy Terms.............................................................. 73

3.3

Generalized Forces due to Spring or Damper ........................................................ 74

3.4

Generalized Forces due to the Control Arm Length Constraint.............................. 75

3.5

Generalized Forces due to the Control Arm Orthogonality Constraint ................... 76

4.1

Kinetic and Potential Energy Terms for the Motion of the Rear Suspension .......... 86

4.2

Kinetic Energy Terms for the Rotation of the Rear Wheels and Tires .................... 87

4.3

Generalized Forces Associated with a Rear Spring or Damper .............................. 87

4.4

Constraint Forces Associated with the Panhard Rod.............................................. 88

7.1

Desired Longitudinal Force Sign and Sign of the Longitudinal Slip...................... 124

9.1

NCSU Legends Car - Front Suspension Geometric Data..................................... 151

9.2

NCSU Legends Car - Rear Suspension Geometric Data, Spring Data and Damper Data ............................................................................ 152

9.3

NCSU Legends Car - Sprung Mass Geometric Data ........................................... 153

9.4

NCSU Legends Car - Model Mass and Inertia Properties .................................... 154

9.5

NCSU Legends Car - Suspension Spring and Damper Properties ........................ 155

viii

9.6

Miscellaneous Tire Model Parameters: Geometric Data, Slip Equation Parameters and Normal Force Characteristics ..................................................... 157

9.7

Delft ’97 Tire Model Parameters: Pure Longitudinal Slip Equation ..................... 158

9.8

Delft ’97 Tire Model Parameters: Pure Lateral Slip Equation .............................. 159

9.9

Delft ’97 Tire Model Parameters: Combined Slip Equations................................ 160

9.10 Driver Model Parameters .................................................................................... 163 9.11 Vehicle Setup Parameters ................................................................................... 164 9.12 Optimization and Cost Function Computation Parameters................................... 165 9.13 Vehicle Suspension Parameter Optimization Ranges ........................................... 177 9.14 Vehicle Suspension Parameter Optimization Results ........................................... 178

ix

List of Figures 2.1

Earth Fixed and Vehicle Sprung Mass Coordinate Systems ................................... 48

3.1

Schematic Showing Front of Sprung Mass and Control Arms ............................... 57

3.2

Schematic of Spindle and Control Arms ................................................................ 58

3.3

Schematic of a Generic Control Arm..................................................................... 70

5.1

Schematic of the Rack and Pinion Steering System ............................................... 90

5.2

Relationship between the P and S Coordinate Systems .......................................... 91

5.3

Relationship between the D and P Coordinate Systems ......................................... 92

5.4

Schematic of the Four Bar Linkage Steering System ............................................. 94

7.1

The Tire Model Coordinate System .................................................................... 118

7.2

Relationship between Tire Velocity Components................................................. 123

8.1

Driver Path for the Kenley, NC Race Track ........................................................ 133

8.2

Steering Profile for the Kenley, NC Race Track .................................................. 135

8.3

Driver Speed Controller Block Diagram.............................................................. 139

8.4

Effect of the Traction Control Gain Parameter on the Acceleration ..................... 143

8.5

The Simplified Preview-Follower Control System ............................................... 146

x

9.1

Schematic of the Kenley, NC Race Track............................................................ 161

9.2

Optimized Steering Profile for the Kenley, NC Simulation................................... 168

9.3

Optimized Velocity Profile for the Kenley, NC Simulation .................................. 169

9.4

Comparison of the Prescribed Velocity and the Actual Vehicle Velocity.............. 170

9.5

Vertical Acceleration of the Sprung Mass (Sprung Mass Coordinate System)...... 171

9.6

Longitudinal Wheel Slip Percentages .................................................................. 172

9.7

Vehicle Position and 9.0 Seconds (Exiting Turn 2).............................................. 173

9.8

Tire Normal Loads.............................................................................................. 174

9.9

Vehicle Lateral Position Error............................................................................. 175

9.10 Yaw Velocity...................................................................................................... 176 9.11 Steering Wheel Angle and Lateral Acceleration................................................... 177

xi

1 Introduction 1.1 Motivation for the Study The demands imposed by the optimal design process form a unique set of criteria for the development of a computational model for vehicle simulation. Due to the large number of simulations which must be performed to obtain an optimized design the model must be computationally efficient. For a fixed execution time a faster simulation will, in general, lead to a better design. A competing criteria is that the computational model must realistically model the vehicle. Current trends in vehicle simulation codes have tackled the problem of realism by constructing elaborate full vehicle models containing dozens if not hundreds of distinct bodies. Each body in a model of this type is associated with six degrees of freedom. Numerous constraint equations are applied to the bodies to represent the physical connections. 1 While the formulation of the equations is not particularly difficult, and in fact has been automated in several software packages, the resulting model requires a considerable amount of computational time to run. This makes the model unsuitable for the application of computational optimal design techniques.

1

Details of this approach can be found in P. Nikravesh’s “Computer-Aided Analysis of Mechanical Systems”.

1

Table 1.1 - Vehicle Model Degrees of Freedom

Body sprung mass rear suspension front right suspension

Degrees of Freedom 7 7 7

front left suspension

7

wheels

4

Constraint Constraint Equations Type 1 Euler Parameter normalization 5 EP norm, panhard rod, trailing links 5 EP norm, upper (2) and lower (2) control arms 5 EP norm, upper (2) and lower (2) control arms none none

Past research in the field of vehicle dynamics has produced numerous computational models which are small enough and fast enough to satisfy the speed demands of the optimal design process. These models typically use less than a dozen degrees of freedom to model the vehicle. They do a good job of predicting the general motion of the vehicle and they are useful as design tools but they lack the required accuracy for optimal design. A model which bridges the gap between these two existing classes of models is required for optimal design. This type of model combines element of both approaches to obtain an accurate solution and yet still emphasize computational efficiency. This is the type of model which is developed in this thesis. The model consists of eight bodies which represent the sprung mass, the rear suspension, the left front spindle, the right front spindle, and the four wheels. There are a total of twenty-eight dynamical degrees of freedom which are distributed as shown in Table 1.1. The total number of computational degrees of freedom is summarized in the Table 1.2. The equations of motion are second order which means that for each dynamical degree of freedom there are two computational degrees of freedom (obtained in

2

Table 1.2 - Computational Degrees of Freedom

Body

sprung mass rear suspension front right suspension front left suspension wheels

Constraint Equations

Dynamical Degrees of Freedom 7 7 7 7 4

1 5 5 5 none TOTAL

Computational Degrees of Freedom 15 19 19 19 8 80

converting the second order differential equations to pairs of first order differential equations). The constraint equations introduce additional degrees of freedom in the form of Lagrange multipliers which are necessary for determining the constraint forces. There are a total of 80 computational degrees of freedom.

1.2 Historical Background 2 The study of automobile stability and control is a relatively new field. Although significant quantities of automobiles were being produced in the early 1900s few efforts were made to quantify the handling issues. Much of the early development was done on a “cut and try” basis and this methodology is reflected in the literature. The majority of the effort prior to 1925 was expended in designing suspensions which would keep the tires in contact with the ground as much as possible in order to enable more effective steering control. This preoccupation with controllability is typical of the early work. Progress in the area of automotive stability was not seen until the 1930s.

2

Much of the historical information prior to the mid 1950s is from the following references: [ Segel, 1956a], [Milliken, 1956], [ Segel, 1956b] and [Whitcomb, 1956].

3

In 1903 the Wright brothers successfully built their first airplane. In the same year G. H. Bryan started his pioneering work on a mathematical theory of airplane stability which was a few years later [Bryan, 1911]. While the refinement of Bryan’s stability theory progressed steadily similar theories for the automobile didn’t appear until much later. This delay was most likely due to the less pressing need to consider stability in ground vehicles as compared to aircraft. The development of usable aircraft hinges on an understanding of aerodynamics and how it affects the stability of an aircraft. This understanding had been evolving with the use of scale models and wind tunnels. The slow development of an automotive stability theory was also the result of a lack of understanding of the role of the tire mechanics in the stability of an automobile. The emphasis on vehicle control between 1900 and 1930 led to kinematic studies of suspension and steering geometries. These studies led to improved designs including Akermann steering geometry. Much of the remaining development work was concerned with the drivetrain, structure and performance of the automobile with one notable exception: A general theory of ride dynamics (the motion of the automobile in its plane of symmetry) was well established by 1925. However, very little, if any, progress had been made in the areas of static and dynamic directional stability. This statement may seem a little strange at first given that the equations for ride dynamics are similar to those involved in a full stability analysis. The key difference lies in the need to understand the mechanism of lateral force generation by the tire. Without this knowledge it is impossible to obtain meaningful stability results.

4

In 1925 Broulheit published a paper in which the basic concepts of side-slip and slip angle were recognized for the first time [Broulheit, 1925]. The recognition of these concepts came about during attempts to explain the phenomenon of steering shimmy which plagued vehicles of the time period. In 1931 Becker, Fromm, and Maruhn published a text on the role of the tire in steering system vibrations and further developed the field of tire mechanics [Becker, 1931]. This realization enabled further study of the problem of automotive stability. During the 1930s the Cadillac Suspension Group of General Motors, under the direction of Maurice Olley, developed the first independent front suspension used on an American car. It was found that certain steering geometries led to a condition which the group termed oversteer [Olley, 1937]. It was recognized that these geometries led to vehicles which were unsafe at high speeds. Olley’s oversteer is recognized today as being roll oversteer. Further investigation revealed that behavior similar to oversteer could be induced by over loading or under inflating the rear tires. In 1934 Olley wrote an unpublished report containing his findings and in which the proposition of oversteer / understeer was stated and the idea of critical speed was first mentioned [Olley, 1934]. As a result of this research Goodyear Tire and Rubber Company began rolling drum tests to determine tire characteristics and in 1935 R. D. Evans published the results in a paper on lateral tire characteristics [Evans, 1935]. This paper gave data on cornering force and selfaligning torque. This work precipitated a period of extensive research at General Motors. The concepts of steady-state directional stability and roll steer were explored. Further

5

exploration of steady-state tire characteristics occurred and skidpad tests were used for the first time. A fundamental understanding of the steady state tire characteristics was developed and a qualitative understanding of the transient behavior was obtained. During the period from 1939 to 1945 very little progress was made due to World War II. In 1950 Lind Walker summarized the current state of knowledge on the issue of directional stability and introduced the concept of the ‘neutral steer line’ and the ‘stability margin’ [Walker, 1950]. These concepts had already been established in aeronautical circles and were suggested as criteria for steady state directional stability in automobiles. The concept of using aerodynamics and tire characteristics to aid in achieving stability was also proposed.

Vehicle Modeling By the middle of the 1950s the groundwork for a mathematical model of the vehicle had been laid. A basic understanding of the tire enabled the creation of reasonably accurate mathematical tire models. In 1956 William F. Milliken, David W. Whitcomb, and Leonard Segel of the Cornell Aeronautical Laboratory, published the first major quantitative and theoretical analysis of vehicle handling in a series of papers [Segel, 1956a][Milliken, 1956][Segel, 1956b][Whitcomb, 1956]. These papers formed the basis for research in the area of automotive stability and control for the next three decades and are still frequently referenced in the current literature.

6

Milliken’s paper [Milliken, 1956] provides a historical overview of the field from which much of the above material was taken. In summarizing the progress made to date, Milliken made the following statement, Thus, [the] major effort in handling research to date has been in the recognition of individual effects, their isolation, and examination as separate entities. This work naturally started out as qualitative and in some instances has become quantitative. It has been conceptual in character; it has been pioneering and not infrequently intuitive and inspired, but it can hardly be viewed as an end in itself. Rather, it is a substantial beginning. All the individual effects now known need quantitative analytical expression. More significant, however, is the need for comprehensive, integrated analysis methods, for such overall theories will enable the prediction of the actual motion by rationally and simultaneously taking into account all of the separate effects.

Milliken also noted that, although a great deal of progress had been made in understanding the tire, there was much to be done still. Although much has been learned about tire modeling Milliken’s observation is still true today; dynamic data on tires is only now becoming available. There were no universally accepted set of reference axes and measured tire data of the period were typically confined to two or three of the possible six force/moments. This made translation of the data from one set of axes to another difficult if not impossible. It was also recognized that the effects of tire design on handling were largely unknown and that there existed a need to perform testing on a wide variety of common passenger car tires to determine the effects of the various design parameters. In discussing of the future objectives of the Cornell Aeronautics Lab research program Milliken emphasized the need to concentrate on the ‘objective analysis of car stability and

7

control’. In the process he made the following distinctions between stability and control, performance and ride. In general, an automobile has ‘six-degrees-of-motion’freedom, and stability and control may be thought of as those lateral motions out of the plane of symmetry involving rolling, yawing and sideslipping. (‘Performance’, by way of distinction, is concerned with fore-and-aft motions in the plane of symmetry, such as acceleration, speed, and braking, while ‘ride’is composed of the vertical and pitching motions in the plane of symmetry.)

The second paper of the series, written by Leonard Segel, derives a set of nondimensionalized linearized three degree of freedom equations for lateral and directional motion [Segel, 1956b]. In accordance with the research goals outlined by Milliken the emphasis of the model was put on modeling for analysis of stability and control. The bounce and pitch degrees of freedom of the chassis were ignored and a fixed longitudinal roll axis parallel to the ground was used. Segel also made several other simplifying assumptions including constant forward velocity, fixed driving thrust divided equally between the rear wheels, and that the lateral mechanical properties of the tires are decoupled from the longitudinal mechanical properties at the speeds studied. The unsprung mass was modeled as a single non-rolling lumped mass. An experimental validation of the model was performed using a 1953 Buick Super fourdoor sedan. The vehicle was put through both pulse steering input and step steering input tests and the transient response for the three degrees of freedom included in the model (lateral displacement, yaw and roll) were measured at a variety of constant forward velocities. The theoretical predictions of the model are compared to the experimental data

8

taken at 32 mph in a series of frequency response curves with the results showing good correlation. The final paper in the series, written by D. W. Whitcomb, draws a series of conclusions on automobile stability and control using a two degree of freedom model (yaw and side-slip) with experimentally determined parameters [Whitcomb, 1956]. Due to the lack of a roll degree of freedom, Whitcomb was able to assume that the car has no width and that the tires lay on the centerline of the vehicle (a “bicycle model”). A set of linearized differential equations is derived using stability derivatives and the steady state and transient responses are studied. In studying the yaw response of the vehicle at a constant vehicle side-slip angle (same angle for both tires) he introduces the concept of the “static margin”. The static margin is an indication of the sense and amplitude of the yawing moment associated with the total tyre side force. It immediately determines the yawing moment that the tyres would provide in reacting an externally applied side force.

In his summary of response characteristics Whitcomb recognizes the strong influence of the static margin on vehicle stability. For vehicle with a negative static margin it was recognized that a critical speed existed, which if exceeded, would lead to instability. As noted by Milliken, there existed a need to quantify and refine the current knowledge of the individual vehicle subsystems. Additionally, he recognized the need to combine these refined models into improved full vehicle models. Progress towards achieving these goals began to be made with research done in the early 1960s.

9

In 1960 H. S. Radt and W. G. Milliken Jr. explored the motions of a skidding automobile [Radt, 1960a][Radt, 1960b]. They used a relatively simple vehicle model with yaw and lateral velocity as the only degrees of freedom. A tire model was incorporated which included the effect of saturation of the side force in the presence of braking and thrust forces via the concept of a friction circle. Results were presented for a series of steady state and transient maneuvers on a low friction surface ( µ=0.3). A simple driver control was also implemented to study skid recovery. The driver model was based on feedback from heading angle with a first order lag. Results are presented for several gains and lag time constants. In August of 1961 Martin Goland and Frederick Jindra published a paper which they used a two degree of freedom (yaw and sideslip) vehicle model to study the directional stability and control of a four wheeled vehicle [Goland, 1961]. The model is a simplified version of Segel’s model with the main difference being that the roll degree of freedom enters as a quasi-coordinate which is only used to calculate the vertical load on the tires. The paper takes into account the effects of load transfer and the variation of the cornering performance of the tires with vertical loading. Results are presented which show how the stability of a vehicle changes as the center of mass is moved, the tire inflation pressure is changed, and the tire tread width is changed. The effect of tread width and inflation pressure on the tire properties is given by a simplified form of the semi-empirical equations published by R.F. Smiley and W. B. Horne in the late 1950s [Smiley, 1958]. Walter Bergman published a paper in 1965 in which he explored the nature of vehicle understeer and oversteer. While the definitions of the terms were relatively well 10

established for steady state maneuvers, they were not well established for the transient case. Bergman discussed the many origins of understeer and oversteer behavior including steering inputs, aerodynamic forces and inertia forces in the transient case. He noted that understeer and oversteer could be recognized by considering the change in the yaw velocity induced by a change in lateral acceleration. This definition is in accordance with the standardized definitions of oversteer and understeer put forth by the Society of Automotive Engineers [S.A.E., 1965]. Bergman also develops a six degree of freedom vehicle model to explore understeer and oversteer behavior as well as vehicle stability. The model consists of a sprung mass and a single unsprung mass. The position of the unsprung mass is given with respect to an inertial coordinate system by a two dimensional vector and a yaw angle. The location of the sprung mass is given relative to the unsprung mass in terms of four vertical wheel displacements. Both masses are assumed rigid which implies that one of the vertical displacements is redundant. In 1966 Segel published a paper in which the stability of a free control automobile (i.e. a vehicle with torque input at the steering wheel as opposed to a steering angle input) was studied [Segel, 1966]. He proposed a two degree of freedom quasi-linear (due to Coulomb friction) model for the steering system. This steering model was added to his three degree of freedom model which was discussed above. The model was validated by comparing simulation output, performed on an analog computer, to experimental data. A reasonably good correlation was demonstrated as long as the lateral acceleration of the vehicle did not exceed 0.3 g. He was able to effectively model the stable and unstable

11

vibrational modes of the combined vehicle and steering model and to relate them to vehicle design parameters. In 1967 R. Thomas Bundorf of General Motors published a paper relating vehicle design parameters to the characteristic speed and to understeer [Bundorf, 1967]. This paper utilized the definitions of understeer and characteristic speed proposed by the SAE publication Vehicle Dynamics Terminology [S.A.E., 1965]. Methods are proposed to predict understeer quality in vehicle designs and for measuring understeer in existing vehicles. It is noted that the characteristic speed is an attribute associated with a linear vehicle model. Bundorf argued that under most normal driving conditions, which he characterized as having lateral accelerations below 1/3 g, a vehicle can be accurately modeled by a linear model. This condition led to the construction of a large diameter skid pad at GM for measuring the characteristic speed; it was not possible to reach high enough vehicle speeds for accurate measurement of the characteristic speed on the existing small diameter pad without exceeding the 1/3 g limit on lateral acceleration. Bundorf derived an expression for predicting the characteristic speed of a vehicle given the design parameters. The vehicle model used in his derivation was a bicycle model with Ackermann (no slip) steering. The paper also contains a discussion, written by A. G. Fonda, of Bundorf’s results with several significant contributions and suggestions. D. H. Weir, C. P. Shortwell, and W. A. Johnson published a paper in 1968 which they explored the role of vehicle dynamics on controllability [Weir, 1968]. Their results were obtained using experimental data and simulation data obtained from a model which combined elements of a nonlinear model developed by H. S. Radt in 1964 [Radt, 1964]

12

and Segel’s earlier models. The model consisted of two unsprung masses representing the front and rear suspension assemblies respectively and a single sprung mass representing the body of the vehicle. The dynamics of the vehicle were described by a linearized set of equations in four degrees of freedom (roll of the sprung mass about a fixed axis, lateral velocity, yaw rate, and axial velocity). The three masses were assumed to posses the same yaw rate, axial velocity and side slip velocity. Provisions were made for a stationary tilted roll axis. In accordance with the inclusion of the axial velocity as a degree of freedom, aerodynamic loads on the vehicle, longitudinal tire forces generated by braking and acceleration and rolling resistance were considered. Dynamic data for a number of automobiles made by U.S. manufacturers is also presented and, as an example, the transfer functions for a typical 1960s sedan were calculated. It was noted that the yaw, lateral velocity and roll modes have undamped natural frequencies of approximately 6 rad/sec at 60 mph. The yaw and lateral velocity modes are highly damped and the roll mode is lightly damped. The roll mode damping ratio was found to be approximately 0.2 to 0.3 and it was found to be largely decoupled from the yaw and lateral velocity modes. Increasing vehicle speed tends to lower the vibrational frequencies and decreases damping which leads to a destabilization of the vehicle. By the early 1970s simulations of vehicle dynamics were becoming more complex and realistic. This was primarily due to advances in computing technology. Prior to the 1970s most simulations were performed on analog computers. These machines were capable of solving the vehicle dynamics problems in real time (since the differential equations were modeled by equivalent electrical component networks) in a cost effective

13

manner. Unfortunately it was very difficult to model nonlinear functions of more than one variable on these machines. Since most tire models are nonlinear functions of more than one variable the accuracy of the simulations was compromised by limitations in the computing equipment. The advent of digital computers allowed researchers to create models containing nonlinear functions. This allowed increased realism in the simulations, however, the slow speed of the digital machines (typically 10 to 100 times slower than real time) meant increased computing costs. In the early 1970s researchers designed simulation codes which ran on hybrid computers which combined digital and analog computing hardware [Murphy, 1970][Tiffany, 1970][Hickner, 1971]. The new computers made it possible to run simulations at real time speeds and at the same time include nonlinearities in the model. A number of papers on computing techniques and on models can be found in the literature. A few of the more significant papers are discussed below. In the early 1970s a vehicle dynamics simulation for a hybrid computer was created by the research staff at the Bendix Corporation Research Laboratories [Tiffany, 1970]. The model was based on the ten degree of freedom model created by R. R. McHenry and N. J. Deleys at the Cornell Aeronautics Laboratory for the Bureau of Public Roads [McHenry, 1968]. The BPR-CAL model was improved by adding four spin coordinates for the wheels and three coordinates for the steering system model. The original BPRCAL model had six coordinates for the sprung mass, one vertical coordinate for each front wheel, and one vertical and one rotational coordinate for the rear axle. The steering system model is based on Segel’s model [Segel, 1966]. At the time of publication the model had been partially validated by comparison of simulation results with the CAL

14

model. The model was upgraded in 1971 to include a dynamically accurate model of a four wheel anti-lock braking system [Hickner, 1971]. In 1973 T. Okada et al described in a paper a seven degree of freedom model for vehicle simulation [Okada, 1973]. The model was used to simulate vehicle handling at the first stage of vehicle design. Five of the degrees of freedom were used to model the vehicle (roll, yaw, pitch, lift and lateral position). The remaining two degrees of freedom were used to model the steering system in a manner similar to that proposed by Segel. The vehicle was assumed to move with constant velocity. A tractive force was applied to maintain constant vehicle speed and compensate for the six components of aerodynamic forces which could be applied to the model. A roll axis which moves vertically in accordance with the wheel travel was included. The effects of roll steer, axle steer, caster, camber, toe-in, and so on were approximated by linear functions based on wheel travel, steer angle etc. The simulation could be run in three different modes: straight-running (with lateral “wind gust” disturbances), stationary circular motion (skid pad), and a slalom mode (to predict critical speed). Gyroscopic effects of the wheels were only included in the straight-running simulation mode where vehicle speeds are high. Steady-state motion was assumed in the skid pad simulation which lead to simplification of the equations of motion and an essentially algebraic system of equations for determining the maximum lateral acceleration. For the slalom course simulation transients of the motion of the vehicle were neglected and constant forward speed was assumed. The path followed was assumed to be periodic with length 2*L where L was the distance between the cones. Galerkin’s method was used to solve for the path. The critical speed was taken to be that

15

speed at which a solution could no longer be obtained. Driver response time limitations were considered as well as vehicle limitations in determining the critical speed. In 1973 Frank H. Speckhart published a paper in which he presented a vehicle model containing fourteen degrees of freedom [Speckhart, 1973]. Six degrees of freedom were assigned to the sprung mass, four degrees of freedom were associated with the suspension movement at the four corners of the vehicle, and four rotational degrees of freedom were assigned to the wheels. He used a Lagrangian approach in deriving his equations. Models were presented for several different suspension configurations. The sprung mass was restricted to pivot about a specified roll axis. It is likely that this is done because the suspension models were relatively simple (two dimensional in the case of the front independent suspension) and did not provide a sufficiently accurate representation of the kinematics involved. As digital computers gradually displaced analog and hybrid machines, primarily as a result of economic concerns, it became necessary to create vehicle dynamics models which were completely digital. The combination of the cost of computer time and the slower solution speed of the digital machines made it desirable to create computationally efficient models. In 1973 Bernard published a paper detailing several time saving methods used in the digital vehicle simulation code created for the Highway Safety Research Institute [Bernard, 1973]. He noted that the important sprung mass motions tended to be in the low frequency range (below 2 Hz) and that the significant wheel hop motions tended to be below 10 Hz. This implied that one should be able to integrate the equations of motion

16

with a relatively large time step (0.005 sec) and obtain accurate results. Unfortunately the cycling of brake torque (as in an anti-skid system) could cause rapidly changing spin derivatives for the wheel degrees of freedom. The relatively high frequency motion required a much smaller time step on the order of 0.0001 second. Bernard proposed an approximate method for dealing with the spin degrees of freedom which allowed the use of the larger time step. This improvement in combination with the use of a specially modified predictor-corrector integration scheme which only updated the wheel-hop derivatives during the corrector phase led to a speed improvement of a factor of five. In early 1976 Frederick Jindra published an interim report for the NHTSA describing a vehicle simulation model being created at John’s Hopkins University Applied Physics Laboratory [Jindra, 1976]. The model was called the Hybrid Computer Vehicle Handling Program (or HVHP) because it was run on a hybrid computer. The HVHP model was derived from a refined version of the Bendix Research Laboratories (BRL) model which was discussed above. The HVHP model was used extensively by Calspan in their study on the influence of tire properties on passenger vehicle handling. The HVHP model contained seventeen degrees of freedom distributed as follows: six for the sprung mass, one for the vertical motion of each front wheel, two for the vertical motion and rotational motion of the rear axle assembly, three for the steering system, and four rotational degrees of freedom for the wheels. The program had an option to use an independent rear suspension model. In this case the model contained two degrees of freedom for the vertical displacements of the rear wheels. The steering system model was a lumped mass model consisting of two degrees of freedom to represent rotation of the

17

front wheels about their steering pivots and one degree of freedom for the translational motion of the connecting steering rod and associated mass elements. Friction and compliance in the steering mechanism were included. The rear unsprung mass was assumed to pivot about a point which was constrained to move along the sprung mass vertical axis. This constraint was an improvement over the traditional fixed pivot point and fixed roll axis. No pivot point was assumed for the independent front suspension; the front wheels were assumed to move vertically with respect to the sprung mass. Due to the difficulty in representing nonlinear functions on the hybrid machine, piecewise linear functions were used to describe the spring force, coulomb friction, damping coefficients, roll stiffness, etc. The camber angle, caster angle, and toe angle were specified as functions of the suspension deflection. Compliance coefficients were used to model the change in camber angle and steer angle due to applied forces and moments at the tire. Radial loading of the tire was computed using a point contact model. In 1977 Kenneth N. Mormon of Ford Motor Company presented a paper [Morman, 1977] containing a detailed three degree of freedom model of the front suspension. The model included the effects of lower control arm bushing compliance along the axis of rotation (but not perpendicular to the axis of rotation) and compliance of the ball joints connecting the tie rod ends to the steering knuckle. The model was derived using a standard Lagrangian approach with constraint equations. A variety of displacement type inputs were applied to the model; the results of the simulation matched experimental results fairly well. In the original model all of the spring, dampers and bushings were assumed to be linear. Improvements could likely be made by replacing the linear elements

18

with appropriate nonlinear relations. It was also assumed that the sprung mass of the vehicle forms inertial coordinate system. In 1981 W. Riley Garrot described an all digital vehicle simulation developed at the University of Michigan [Garrot, 1981]. The model contained a total of seventeen degrees of freedom distributed in a manner identical to the HVHP model discussed above. To reduce computational costs the steering system was described statically and the wheel-spin degrees of freedom were handled algebraically. The model contained numerous features which could be turned on or off as desired. These features included an anti-lock braking system, multiple tire models, optional activation of nonlinear kinematic terms, solid rear axle or independent rear suspension and interactive capability. The program was constructed in a modular fashion to enable future enhancements and upgrades. The simulation consisted of two main parts: a vehicle model called IDSFC and a generalpurpose driver module called DRIVER. The driver module could be readily altered without affecting the vehicle model. The driver model controlled steering, braking and drive torque inputs to the vehicle model. It contained five preprogrammed open-loop maneuvers and could accept user defined maneuvers using tabulated data or a user defined subroutine. Various closed-loop control strategies were implemented including a crossover model for path following and two types of preview-predictor models. Mixed open-loop and closed-loop control could be used. Validation of the model was performed by comparison with the validated HVHP model. In 1986 R. Wade Allen and several associates from Systems Technology Inc. performed experimental tests and correlated the results with a computer model in order to

19

validate a simplified lateral vehicle dynamics model and the associated tire modeling procedure [Allen, 1986]. The tests consisted of a number of steady state skid pad runs and several low amplitude sinusoidal steer frequency sweeps while negotiating a steady turn. The tests were performed for a rear wheel driver 1980 Datsun 210 and a front wheel drive 1984 Honda Accord. Several types of tires were used on the Datsun including both radial and bias ply tires. The physical parameters describing the vehicles and the tires were input into a simplified lateral handling model which was derived directly from Segel’s original model [Segel, 1956b] and which was discussed above in the review of D. H. Weir’s paper [Weir, 1968]. A good correlation was obtained with the experimentally obtained data. The model was also used to extrapolate vehicle behavior under combined cornering and braking. In 1987 Allen published a revised model containing five degrees of freedom [Allen, 1987a]. The new model added pitch and forward velocity degrees of freedom and was called VDANL (Vehicle Dynamics Analysis : Non-Linear). It was also a nonlinear model and, unlike earlier linear models, the solutions were obtained in the time domain using numerical integration. Neither of the models approach the complexity of some of the other more detailed models discussed above [Okada, 1973][Speckhart, 1973][Jindra, 1976]; the intent was to provide a simulation code which could be run on relatively inexpensive desktop PCs and which could utilize the graphical output capabilities of those PCs. In 1987 Andrez Nalecz presented the results of an investigation into the effects of suspension design on the stability of vehicles and, in particular, how the design of the suspension related to movement of the roll axis [Nalecz, 1987]. Twenty-five different

20

suspension types were considered. A typical three degree of freedom lateral dynamics model was used with the addition of a quasi-static pitch degree of freedom. The sprung mass was assumed to rotate about a roll axis whose position varied as a function of body roll. The location of the front and rear roll centers was found via a kinematic analysis of the suspension in which the wheel contact patches were treated as revolute joints and were allowed to move laterally along the ground (thus allowing for track width changes). It was found that for certain types of suspensions, most notably the double wishbone and MacPherson strut type systems, that the assumption of a fixed roll axis could not be justified. In 1992 Nalecz published a second paper in which he described an eight degree of freedom model called LVDS (Light Vehicle Dynamics Simulation) [Nalecz, 1992]. The model consisted of a three degree of freedom lateral dynamics model coupled to a five degree of freedom planar rollover model. The models are coupled through the inertia terms and tire force terms. The lateral dynamics model was derived in the same manner as Segel’s original model. The rollover model consisted of sprung and unsprung masses connected through the various elements of the suspension system. The model also included aerodynamic effects; all six possible forces and moments are modeled. The effects of lateral and longitudinal weight transfer were accounted for in determining the lateral forces generated by the tires. The roll axis was modeled in the quasi-static fashion discussed above. In the early 1990s R. Wade Allen and his associates at Systems Technology Inc. published a number of papers in which they validated their VDANL simulation code [Allen, 1992] and in which details of experimental studies and simulation runs involving

21

of vehicle stability and vehicle rollover are presented [Allen, 1990][Allen, 1991][Allen, 1993]. VDANL and IDSFC (which is derived from the HVOSM simulation code) were also put through a rigorous validation process by Gary J. Heydinger et al at Ohio State University [Heydinger, 1990]. Both validations were carried out by comparing experimental data to simulation data in the time domain and in the frequency domain. The control inputs from the experimental tests were recorded along with the vehicle responses for later use as simulation inputs. Sinusoidal frequency sweeps and step inputs were used in the testing. Heydinger explored the use of pulse inputs which require shorter test runs and could excite the same frequency range in a later paper [Heydinger, 1993]. In studying vehicle stability and rollover stability the authors gathered model parameter data for a total of 41 different vehicles of various types. The connection between load transfer distribution (which is largely governed by the relative roll stiffness at the front and rear axles) and vehicle stability was discussed in detail. Simulation results for a set of maneuvers were plotted. The effects of braking, acceleration and throttle lift on stability in limit handling situations was also discussed. A similar paper, also using the VDANL software, was written by Clover and Bernard at Iowa State [Clover, 1993]. Details of the updated vehicle dynamics model VDANL were presented in [Allen, 1991]. The biggest change in the model was the removal of the fixed roll axis assumption and the addition of a front suspension model which reflects camber change with body roll. The model also included the effect of lateral deflection of the tire, wheel and suspension which decreases the track width and affects rollover stability. In [Allen, 1993] the authors demonstrated that the standard single lane change maneuver was sometimes inadequate for vehicle stability

22

studies in that it failed to cause unstable behavior and that it did not adequately model the large lateral displacements which could occur in real world accident avoidance maneuvers. Simulation results for larger lateral displacements (with the same peak lateral acceleration) demonstrated both spinout and rollover. By the early 1980s a shift in the vehicle modeling process was taking place. The demand for accurate vehicle dynamics models combined with the difficulty in deriving the equations of motion for large multibody systems led to the use of general multibody simulation codes. A wide range of capabilities are present in modern MBS codes including the ability to handle non-inertial reference frames, to incorporate flexible elements in the model, to utilize generalized coordinates, and to symbolically generate the equations of motion. Several reviews of multibody codes have been published in recent years, several of which are discussed in more detail below. Additionally, brief descriptions of a few papers utilizing MBS codes for vehicle dynamics simulations are presented below. In 1985 W. Kortüm and W. Schiehlen presented a paper [Kortüm, 1985] which they presented the desirable qualities of an MBS program, discussed two contemporary examples in some detail and utilized the two codes to generate some simple vehicle models. The first code discussed was NEWEUL which generates the equations of motion in symbolic form with the output being FORTRAN code. It had the capability of using both Cartesian and generalized coordinates, non-holonomic constraints and moving reference frames. The second program was MEDYNA which generates the equations of motion in numerical form. It also had the capability of using generalized coordinates and

23

moving reference frames. Both codes supported the use of closed loops (i.e. four bar linkages). In 1993 W. Kortüm and R. S. Sharp published a supplement to the periodical Vehicle System Dynamics in which the capabilities of 27 currently available multibody simulation codes and general purpose vehicle simulation codes were reviewed [Kortüm, 1993]. The programs discussed include ADAMS, MEDYNA, NEWEUL, DADS, AUTOSIM, and SIMPACK among others. Tables were presented which offer comparisons of the capabilities of the various codes. Kortüm discussed the desirable traits of a multibody code and gave a brief discussion of the contemporary numerical methods which are most applicable to vehicle dynamics simulation. Sharp discussed the four models which were used in benchmarking and evaluating the codes in his introduction. In 1994 R. S. Sharp wrote a paper in which he compared the capabilities of the major multibody computer codes with emphasis on those which generate the equations of motion symbolically [Sharp, 1994]. The codes reviewed were selected based on their applicability to automotive simulation. He discussed the methods used by each code in deriving the equations of motion with attention to the limitations of each method. In particular he noted the limitations of each code with respect to the types of constraint equations that could be handled. References to significant papers in the area of multibody dynamics were given. R. J. Antoun discussed a vehicle dynamic handling computer simulation created using the multibody code ADAMS ( Automatic Dynamic Analysis of Mechanical Systems) in a paper which was published in 1986 [Antoun, 1986]. A model of a 1985 Ford Ranger 24

pickup truck was created utilizing a combination of the standard ADAMS model definition language and user written subroutines for non-standard system components such as the tires. A detailed kinematic model of the front I-beam suspension and the rear leaf spring suspension (using a three link approximation) was constructed. The effects of bushing compliance were included in the model. Nonlinear shock absorbers were used. Excellent agreement of simulation results with experimental data was obtained. Other studies were made using models for a 1986 Bronco II and a 1986 Aerostar van. Using the respective models the researchers were able to optimize the stabilizer bar dimensions and tire characteristics at an early stage of the design process. The Bronco model contained 55 degrees of freedom. It was noted that the extensive graphical display capabilities of the ADAMS program were invaluable in debugging the model geometry and in interpreting the results. A paper describing a model built utilizing a program which automates the generation of the equations of motion was presented in 1991 by C. W. Mousseau [Mousseau, 1991]. The program, AUTOSIM, was used to create a 14 DOF vehicle model. The program used a form of Kane’s equations to derive the equations of motion and applies extensive algebraic and programming optimizations to achieve high efficiency. The user was responsible for choosing the generalized coordinates which describe the configuration of the system. It was not necessary to use Cartesian coordinates and numerous constraint equations to formulate the equations of motion. In generating the vehicle model the location and orientation of the spindle was expressed in terms of four generalized coordinates; the generalized coordinates were specified as prescribed cubic

25

polynomial functions of the suspension deflection (a quasi-static approximation). The cubic polynomials were obtained from a kinematic suspension model. The effects of suspension geometry and suspension bushing compliance were included in the suspension model which was also created using AUTOSIM. Integration of the resulting FORTRAN model produced good correlation with measured data. The computational efficiency of the resulting model allowed it to be used in real time in a driving simulator. In 1993 Michael W. Sayers published a paper in which AUTOSIM was used to generate a number of vehicle models [Sayers, 1993]. The simplest model possessed 4 degrees of freedom system while the more complicated models contained 10 degrees of freedom. The emphasis in the paper was on demonstrating the ease with which computationally efficient models can be generated and tested. Yoshinori Mori et al at Toyota described a model created for simulation of active suspension control systems in a paper presented in 1991 [Mori, 1991]. The vehicle model was described using a simulation language. The control algorithms were coded in FORTRAN and interfaced to the vehicle model. The vehicle model contained 20 degrees of freedom. The unsprung masses were assigned three degrees of freedom each and the sprung mass was given six degrees of freedom. Each of the front wheels was assigned a single steer degree of freedom. The model also included a 19 degree of freedom drivetrain model. Provisions for front wheel drive, rear wheel drive and four wheel drive were made. The road surface was modeled using a combination of a flat or undulating surface and random input noise.

26

In 1989 a research group at the University of Missouri-Columbia began a DOT sponsored project to study the effects of vehicle design on rollover propensity [Nalecz, 1988]. A nonlinear 14 degree of freedom vehicle model called the Advanced Dynamics Vehicle Simulation (ADVS) was developed to carry out this research. The model was derived using a Lagrangian approach and utilizes quasi-velocities to describe the angular velocities. The degrees of freedom were utilized as follows: three translational and three rotational for the sprung mass, two for the front suspension and two for the rear suspension and one rotational for each wheel. To study vehicle-terrain interaction it was necessary to model the body of the vehicle as well as the terrain [Lu, 1993]. The vehicle body was represented by a set of massless, three-dimensional nodes which obey nonlinear force-deflection curves. Each node was checked for interference with the terrain at each time step of integration and its position was adjusted as necessary. The force resulting from body-terrain interaction was applied to the vehicle dynamics model. The terrain was modeled by a single curve which was extruded along the direction of travel. This prevented the use of curved roadways and other such fully three-dimensional structures but it simplified the body node-terrain interference calculation substantially. In 1993 the results of a program at Lotus Engineering to develop a vehicle simulation code for studying the application of predominantly linear control algorithms to the suspension of a nonlinear vehicle were published by J. G. Dickinson and A. J. Yardley [Dickison, 1993]. Although commercial multibody simulation codes were available it was desired to utilize a simpler model which did not require the large quantities of descriptive data associated with the more complicated codes. The model which was presented in the

27

paper utilized six degrees of freedom for the sprung mass. The front and rear suspensions were modeled in a quasi-static fashion. Each wheel was assigned a ‘bump’ degree of freedom which was measured relative to the sprung mass. The location of the instantaneous pivot axis was determined from a look-up table based on the value of the bump variable. Since the motion was handled in a quasi-static fashion the pivot axis location, camber angles, wheel hub location, toe angle, effective spring rate, effective damper velocities and so on could be calculated off-line. The front suspension was modeled in the same fashion but adds a steering swivel axis and two degree of freedom steering system. The tires were modeled using the Pacejka curve fits to measured tire data. The longitudinal force at the tires was set by the driver acceleration input. The lateral force was reduced accordingly by utilizing a standard friction ellipse. Wheel angular velocities were apparently not included as degrees of freedom in the model. The authors claimed a speed advantage of a factor of three over more complicated models generated using standard multibody codes and hoped to increase the advantage to a factor of six in later versions of the software. In 1996 Michael R. Petersen and John M. Starkey described a relatively detailed straight line acceleration vehicle model for predicting vehicle performance [Petersen, 1996]. The model included longitudinal weight transfer effects, tire slip, aerodynamic drag, aerodynamic lift, transmission and driveline losses and rotational inertias of the wheels, engine and driveline components. A manual transmission was assumed with 100% clutch engagement. Shifts were simulated by disengaging the clutch completely, assuming that the engine torque is zero during the shift, changing the gear ratio, and then reapplying

28

the full torque of the motor. Shifts occurred when the applied torque at the rear wheels in the next gear exceeded the torque at the rear wheels in the current gear, or alternatively, when redline was reached. After validating the model the authors conducted sensitivity analyses to determine which design parameters most strongly affected vehicle performance.

Driver Modeling Beginning in the early 1960s an increasing emphasis on vehicle safety created a push toward modeling vehicles under the more demanding conditions associated with crash avoidance maneuvers. In order to accurately represent the reactions of the vehicle under these circumstances it was necessary to include the driver as an integral part of the model. While this fact had been recognized in the early 1960s it was not until the late 1960s that increasing computational power and an improved understanding of vehicle dynamics and driver behavior made it practical to model the driver and vehicle together. In 1968 David H. Weir and Duane T. McRuer of Systems Technology Inc. published the first in a long series of papers on modeling driver steering control (lateral control) [Weir, 1968b]. The vehicle dynamics were modeled using Segel’s equations [Segel, 1956b]. The equations were Laplace-transformed and the analysis was performed in the frequency domain. The transfer functions relating the motion variables to the inputs were taken from Weir’s earlier paper [Weir, 1968a]. Although Segel’s steering system model was available [Segel, 1966], the lack of dynamic data on vehicle steering systems made its use impractical. Consequently, a pure gain was used to describe the steering

29

system dynamics. The driver model was divided into four subsystems: quasi-linear compensatory control, pursuit control, precognitive control and a remnant. The quasi-linear compensatory control consisted of a describing function with parameters which were adjusted to fit the situation and the system, an additive remnant and a set of adjustment rules. The form of the driver model, of the describing function and of the parameter adjustment rules was derived from extensive experiments involving human operators. It was noted that the parameter adjustment rules could be eliminated by considering the combined response of the vehicle/driver system. In this case an approximate crossover model was found to represent driver/vehicle behavior adequately. This simplification was a result of experimental studies involving human drivers which found that drivers adjust their behavior to obtain an approximately invariant form for the combined vehicle-driver response function. The pursuit control subsystem modeled the driver’s ability to see the roadway ahead. This is in contrast to the compensatory subsystem in which the driver reacted to errors in the current position of the vehicle. The details of pursuit control are not mentioned except to note that experimental evidence indicates that the magnitude of the feedforward describing function was approximately equal to the inverse of the magnitude of the vehicle response function. Thus the command path and the actual vehicle path are nearly identical. It was also noted that compensatory control was often used in combination with pursuit control to regulate errors in path following. The precognitive control model attempted to mimic learned driver responses. A common example of this type of maneuver is pulling out and pulling in while passing

30

another automobile. Weir and McRuer note that these types of maneuvers do not involve a feedback based on position information or a feedforward based on the desired path. The maneuver is initiated by the driver in response to stimuli other than those involved in pursuit and compensatory control. No other results are presented by Weir and McRuer beyond defining the nature of precognitive control. The driver remnant component of the model accounted for the portion of the driver’s output which was not linearly correlated with the input. It was modeled as a random input which was described by an experimentally obtained power spectral density. It was noted that the major source of this remnant is due to variation of the parameters of the driver describing function. The remnant could be neglected for vehicles which demonstrated good response characteristics. Following discussion of the various model components them authors presented the results of a guidance and control analysis of the potential loop closures for compensatory control. A number of multiloop structures were considered. The best multiloop feedback structures were considered to be those which demonstrated good frequency response and required minimal driver attention. Based on the author’s analysis it was concluded that a feedback structure based on heading angle and lateral acceleration gave the best results. A review of perceptual experiments performed by other authors [Gordon, 1966a][Gordon, 1966b][Crossman, 1966] corroborates Weir and McRuer’s conclusions. In a later paper Weir and McRuer reviewed data from experiments on the directional response of vehicles subjected to cross wind gust disturbances. Driver/vehicle describing functions were measured for several test drivers. The results support Weir and

31

McRuer’s earlier assertion that the driver’s steering outputs could be explained as functions of lateral position and heading angle or alternatively as functions of path angle and path rate. In 1975 Errol R. Hoffman presented a paper [Hoffman, 1975] in which he reviewed the state of the knowledge of human control of road vehicles. He covered lateral control and longitudinal control of automobiles and motorcycles. The areas of research reviewed in the paper were divided into the following major categories: lateral control of automobiles, lateral control of motorcycles, longitudinal control of automobiles and combined lateral and longitudinal control. The relevant portions of Hoffman’s review of the literature in the areas of lateral and longitudinal control of automobiles is summarized below. The work done in the area of lateral control was divided it into four sub-groupings: lateral control vehicle dynamics, perceptual studies, mathematical models of driver steering control and vehicle characteristics and driver/vehicle performance. Hoffman classified the work done on lateral control vehicle dynamics category into the following three categories: fixed control, free control and vehicle-driver interface variables. Fixed control occurs when steering wheel input angle is specified directly. Free control occurs when the steering wheel input is in the form of a specified torque. The majority of the research up to the time of publication of Hoffman’s review had been performed on the fixed control mode; very little work had been performed using free control. Hoffman noted that, in reality, a human driver uses a combination of the these two types of control. He also noted that the proportion of each type of control varies with the type of maneuver being performed. Hoffman’s third grouping under lateral control

32

vehicle dynamics category encompassed research done on driver/vehicle interface variables. Driver/vehicle interface variables are defined as the quantities which relate steering wheel input (either angle or torque) to vehicle response. Typically they are approximations to the actual output and are used in determining gains in the control algorithms. Again, the majority of the existing work concentrated on identification of the gains associated with fixed control (i.e. neutral steer path curvature vs. steering wheel angle, etc.). Very little work had been performed relating steering force to vehicle response for the free control mode. At the time of Hoffman’s review a number of papers on driver perception of the roadway had been published. Several papers suggested that the driver uses the perceived velocity field to guide the vehicle. Later studies of driver eye movements indicated that peripheral vision is used to monitor steering control for tracking and directional guidance while central vision is used for obstacle avoidance. Studies of driver steering control movements indicated that vehicle yaw rate and inertial lateral deviation are the most probable control cues used by the driver. Hoffman reviewed a variety of mathematical models available in the literature at the time. He included a brief review of the quasi-linear model proposed by McRuer et al which was discussed in above. He also briefly discussed the predictive models of Kondo and Ajimine and of Yoshimoto [Kondo, 1968][Yoshimoto, 1969]. These models were single loop models which used estimated position and heading data as feedbacks. Hoffman also mentioned an optimal control model outlined by Roland and Sheridan [Roland, 1966][Roland, 1967] which was useful in course planning situations. He noted that the

33

primary difficulty in using this type of model lies in determining the driver’s cost weighting. Several other types of models were briefly reviewed. The research done in the area of vehicle characteristics and driver/vehicle performance was primarily concerned with relating driver/vehicle response to vehicle design parameters and vice versa. A considerable amount of work had been performed at the time of Hoffman’s review. The majority of the papers were concerned with determining the optimal vehicle characteristics which maximized driver/vehicle dynamic performance. Steering force, stability factor, and steering gain and sensitivity were among the design parameters considered. Hoffman noted that the state of the knowledge for longitudinal control was considerably less developed and that very little work had been performed for combined lateral and longitudinal driver control models. A good deal of work had been presented regarding the drivers perceptual processes for longitudinal control and a number of driver models had been created. Some work related to the car-following control task had been presented. Combined models containing both lateral and longitudinal driver controls were essentially nonexistent. The lack of full vehicle models which included the interface variables for braking and acceleration largely prevented the application of driver models to vehicle dynamics simulations. In 1977 Duane T. McRuer et al presented a paper which reviewed the progress made to date in driver modeling with emphasis on quasi-linear models for lateral control [McRuer, 1977]. He noted that the driver can be modeled in terms of three types of controllers which, based on the position and velocity of the vehicle and information from

34

the road ahead of the vehicle, regulate steering wheel angle. The three controllers were the same ones introduced in earlier papers by the same authors: compensatory, pursuit and precognitive. Under normal driving circumstances any one of these control modes could be used and sometimes they were combined. In the later part of the paper McRuer presented some results which verified some of the assumptions made in the paper and assisted in quantifying the model. In the same issue George A. Bekey presented a paper [Bekey, 1977] in which he discussed a variety of driver models used for car following tasks. The models assumed that the positions and velocities of the lead car and of the following car were known by the driver. The first part of the paper was devoted to discussion of models using classical control structures. It was found that, in general, linear models performed quite well for small disturbances in the neighborhood of the steady state condition. The nonlinear models functioned better than the linear models and they also performed well under transient conditions. The second part of the paper discussed models derived from optimal control approaches to the problem. A quadratic criterion function was used which was based on spacing error, velocity error and control effort. The control algorithm was determined by minimizing this function. The resulting control model did not include driver reaction time, neuromuscular dynamics or vehicle nonlinearities. Incorporation of time delays and vehicle nonlinearities into the model and experimentally determining the gains provided improved the results. Optimal stochastic models which included the effects of input and output noise provided additional improvement. Two other classes of models were reviewed: look ahead models and finite state models. The look ahead model was similar to the above models but

35

with the addition that the driver “averages” the behavior of several lead vehicles to determine a course of action. The finite state model arranges the following task into four specific conditions. Driver action was based on the current condition and a set of rules was formulated for switching modes. Both of these models provided acceptable results. In 1978 David H. Weir and Richard J. DiMarco presented a paper [Weir, 1978] in which the results of vehicle handling tests from three sources were correlated to vehicle design parameters and then evaluated quantitatively and subjectively. The original tests were conducted by Systems Technology, Inc. (STI), the Highway Safety Research Institute (HSRI), and the Texas Transportation Institute (TTI). The STI study included data obtained with an expert test driver as well as data from typical drivers. The tests included unexpected obstacle avoidance maneuvers, lane change maneuvers and double lane change maneuvers. The results were correlated with a simplified two degree of freedom analytical model in order to obtain estimates of the vehicle parameters. The drivers also gave subjective evaluations of the handling of the various vehicles and these results were used to determine the limits of acceptability on the vehicle design parameters. The results for the expert test driver are analyzed and presented separately. Edmund Donges wrote a paper in 1978 in which he described a two level model of driver steering behavior [Donges, 1978]. The model consisted of a compensatory submodel in parallel with an anticipatory submodel (which was equivalent to McRuer’s pursuit submodel). The mathematical forms of the submodels were derived using parameter identification techniques. The experimental data was obtained using a driving simulator and a number of typical test drivers. The vehicle dynamics model consisted of

36

two components. The first component modeled the lateral dynamics of the vehicle and the second modeled the longitudinal dynamics of the vehicle. Both components were extremely simplified. The parameter estimates of both of the driver submodels showed significant dependence on vehicle speed and on the curvature of the roadway. Comparison of the performance of the simulated driver and of the real drivers showed good agreement. In 1979 R. Wade Allen and Duane T. McRuer expanded on their compensatoryprecognitive-pursuit quasi-linear model [Allen, 1979]. Previous papers had concentrated on the compensatory and precognitive aspects of the human controller. This paper analyzed more recent driving simulation results from a pursuit control point of view and proposed a pursuit control model with roadway preview. The curvature of the roadway was provided as an input in addition to the standard heading angle and lateral position inputs. During the 1980s R. Wade Allen, Henry T. Szostak, Theodore J. Rosenthal, et al published a series of papers [Allen, 1982, 1986, 1987a, 1987b] based on work done for the NHTSA [Allen, 1988]. The first paper, published in 1982, reviewed and expanded on previously published driver steering control models. Experimental data were used to validate the model structure. The second paper, published in 1986, presented test methods and modeling procedures for identifying the directional handling characteristics of vehicles. Computer modeling of the vehicle dynamics was used to extrapolate vehicle response beyond the typical steady-state tests done in previous papers. The third paper, published in 1987, presented linear dynamic models, nonlinear dynamic models and numerical procedures to implement them on a typical microcomputer. Both front and rear

37

wheel drive vehicles were analyzed in the paper. The fourth paper, also published in 1987, presented an updated driver control model. The primary feedback was based on perceived curvature error in the vehicle path with a secondary feedback to control lateral position error. Simulation results obtained by combining the updated driver model with a modified version of the nonlinear vehicle model from the previous paper are presented. The simulation results from several accident avoidance type maneuvers were presented including results which demonstrated the transition to oversteer under combined braking and cornering. In September of 1993 A. Modjtahedzadeh and R. A. Hess published a paper which presented a simple model of driver steering control [Modjtahedzadeh, 1993]. The model was capable of producing driver/vehicle steering responses which compared favorably with experimental data. The results of a computer simulation of a lane keeping task on a curved roadway were provided for both two-wheel steer vehicles and four-wheel steer vehicles. In 1993 K. Guo and H. Guan published a brief review and quantitative comparison of the various approaches used in implementing driver/vehicle/road closed-loop direction control systems [Guo, 1993]. The review includes a discussion of compensation tracking models and preview tracking models. Utilizing optimal control theory, the authors develop a theoretical framework for designing an optimal preview controller for use with vehicle simulation codes. The resulting controllers are compared to the other types of controllers reviewed in the paper.

38

In 1996 R. Wade Allen et al of Systems Technology Inc. published a paper which discussed an updated driver control model [Allen, 1996]. The lateral position control model was updated to include a feedforward based on curvature error. This addition to the model enhanced stability at higher lateral accelerations. The paper also discussed a driver speed control model. When the roadway is straight the driver model follows a preprogrammed speed profile (i.e. speed limits). In the event that a curve lies ahead the speed control model is capable of utilizing the road curvature information to maintain safe or comfortable speeds through curves, reducing speed if necessary, based on a user selected maximum lateral acceleration value. Model responses were shown for several test cases which demonstrated improved controller stability over earlier models. The lateral position control model exhibited some oscillatory instability as the vehicle approached its cornering limits.

Model Parameter Measurement As vehicle models have increased in complexity it has become more difficult to obtain accurate values for the parameters which describe the model including inertial parameters and suspension geometry parameters. Additionally, for the more complex vehicle models, it is necessary to obtain sufficient information to describe environmental inputs to the model such as road surface profile and wind velocity information. There are two approaches to determining these parameters: direct measurement and parameter identification. Modern simulation codes should use a combination of these methods. Several research papers on these approaches are described below.

39

In addition to the difficulties associated with measuring parameters there is the problem of validating the input. The large number of parameters which describe even a moderately complex vehicle model can make it difficult to discover errors. Recently two papers were published which discuss this problem [Bernard, 1994][Gruening, 1996]. The authors note that there are several potential sources of error including erroneous measurements, misinterpretation of the input parameter format, and typographical mistakes in data entry. To detect these errors it is suggested that simulations be run for a series of standard maneuvers and that the results be checked against closed form solutions associated with simpler models (this only applies to low lateral acceleration maneuvers which exhibit essentially linear behavior). Additionally, the authors suggest that a number of basic vehicle performance metrics be computed. While these methods may not be completely foolproof they can reliably detect errors in the input data and in the model itself. Direct Measurement A number of papers have been published dealing with measuring various vehicle model parameters. Several road profile measurement systems have been developed which allow data to be acquired while traveling at freeway speeds. Tetsushi Mimuro et al describe a system of this type which use four laser displacement transducers attached to the bottom of a vehicle [Mimuro, 1993]. Several devices for measuring suspension parameters have been developed to facilitate vehicle modeling. Some require removal of the suspension components from the vehicle (e.g. [Bell, 1987]) while others can obtain measurements without disassembly as done in [Chrstos, 1991]. A facility for the

40

measurement of vehicle inertial parameters is described in [Heydinger, 1995]. The system is capable of measuring center of gravity location, roll, pitch and yaw inertias and the rollyaw product of inertia. Parameter Identification The application of parameter identification techniques to vehicle models is a recent development. A search of the literature revealed a total of four papers on the subject. A literature search performed by the authors of one of the papers found even fewer papers on the subject. Brief reviews of the papers are presented below. Richter, Oberdieck and Zimmerman appear to have published the first paper in which a form of parameter identification was used to improve the fit of a vehicle model to experimental data [Oberdieck, 1979]. In their paper they utilized a simple bicycle model with a nonlinear tire side force coefficient model in order to improve model performance. All of the model parameters are easily obtained except for the side force coefficients which include axle kinematics and elasticities in addition to the tire characteristics. The nonlinear side force model contained three unknown constants for each axle and a seventh parameter which specified the maximum obtainable lateral acceleration. A least squares fit method was used in combination with a genetic search algorithm to obtain the best fit to the experimental data. For the purpose of demonstrating the technique the authors utilized a complex 27 degree of freedom model to generate the ‘experimental’ data instead of instrumenting a real vehicle. Y. Lin and W. Kortüm published a paper in 1991 on parameter identification in the time domain [Lin, 1991]. They presented a method which is applicable to linearized

41

models containing possibly nonlinear forcing terms. It was assumed that the model is linear in the unknown parameters. A least-squares type cost function was utilized and a closed form solution was obtained for the optimal set of parameters. The authors presented an example of the technique using a four degree of freedom bicycle model. A simulation was run using a set of parameters selected by the authors and a colored noise input. The parameter identification algorithm was then applied to the model using the specified input and the simulation output. The parameter identification yielded results which correlated very well with the known values. In 1993 Feng Huang et al. published a paper in which parameter estimation was used to determine model parameters [Huang, 1993]. The authors of this paper were only able to locate one other author who had previously published a paper on using parameter estimation in vehicle handling dynamics [Oberdieck, 1979]. Experimentally obtained data from sinusoidal sweep steering tests was reduced to obtain the frequency response functions between the various input and output variables. The parameter estimation process involved fitting the model responses in the frequency domain to the experimentally obtained frequency response functions. Both models used in this paper were linear which restricted their use to situations in which the lateral acceleration was small. The first model was a simple two degree of freedom bicycle model. The second model added an artificial roll axis and tire camber effects and contained three degrees of freedom. Yi Kyongsu and Karl Hedrick presented a paper in mid-1993 in which they used a sliding observer to estimate unmeasured states during the parameter identification process [Kyongsu, 1995]. The authors consider nonlinear systems of equations which are linear in

42

the unknown parameters. A sufficient condition for convergence of the algorithm was given. In the example given the observer was derived from a quarter car model and was used to identify the parameters of a half car (bicycle) model. Experimental data was obtained from a laboratory based half-car test rig and was used to validate the parameter identification process. Comparison of the measured damper force-velocity relationship and the identified force-velocity function shows good agreement.

1.3 Derivation Methodology and Overview of the Thesis The equations of motion of the vehicle are derived using a slight variation on Lagrange’s formulation. Lagrange’s equation is typically written

d  ∂L  ∂L  − = Qk dt  ∂q&k  ∂qk

(1.3.1)

L ≡T − V where L is Lagrangian and Qk is the generalized force associated with the generalized coordinate qk. The Lagrangian L is defined as the difference between the kinetic energy of the system (T) and the potential energy of the system ( V). The generalized forces are those forces not included in the equation via the potential energy term V due to their nonconservative nature. These forces consist of the forces generated by the tires, by the aerodynamics of the vehicle, and by the suspension dampers. The potential energy of the system is the result of gravitational potential energy and the energy stored in the suspension springs and tires (modeled as springs in the vertical direction). Substituting the definition of the Lagrangian into the equation gives the result

43

d  ∂T  d  ∂V  ∂T ∂V + = Qk  −  − dt  ∂q&k  dt  ∂q&k  ∂qk ∂qk

(1.3.2)

The potential energy associated with the vehicle depends only on the position of the bodies (i.e. the qk) which implies that the second term of the equation above must be zero. d  ∂T  ∂T ∂V + = Qk  − dt  ∂q&k  ∂qk ∂qk

(1.3.3)

The equation can be broken down further by dividing the kinetic and potential energy terms into sub-terms associated with the various subsystems of the vehicle model. T = TSM + TRS + TRT + TLF + TRF + TFT V = VSM + VRS + VLF + VRF Qk = QRS,k + QLF,k + QRF,k + QT,k + QC ,k + QRP,k + QTR ,k + QCA ,k

(1.3.4)

The sub-terms are identified in Table 1.3. Substituting these expressions into Equation 1.3.2 and rearranging gives d ∂TSM  d  ∂TRS  d  ∂TRT  d  ∂TFL  d  ∂TFR  d ∂TFT   +  +  +  +  +   dt  ∂q&k  dt  ∂q&k  dt  ∂q&k  dt  ∂q&k  dt  ∂q&k  dt  ∂q&k  −

∂TSM ∂TRS ∂TRT ∂TFL ∂TFR ∂TFT ∂VSM ∂VRS ∂VFL ∂VFR − − − − − + + + + ∂qk ∂qk ∂qk ∂qk ∂qk ∂qk ∂qk ∂qk ∂qk ∂qk

− QRS,k − QLF,k − QRF,k − QT,k − QC,k − QRP,k − QTR ,k − QCA ,k = 0

44

(1.3.5)

Table 1.3 - Identification of Sub-Terms in the Equations of Motion

Term TSM TRS TRT TLF TRF TFT VSM VRS VLF VRF QRS,k QLF,k QRF,k QT,k QC,k QRS,k QTR,k QCA,k

Description Kinetic energy of the sprung mass Kinetic energy of the rear suspension including the terms associated with the mass of the rear wheels and tires. Kinetic energy associated with the rotational motion of the rear tire & wheel assemblies Kinetic energy of the left front suspension Kinetic energy of the right front suspension Kinetic energy associated with the front wheels and tires. Gravitational potential energy of the sprung mass Gravitational potential energy of the rear suspension including the rear tires and wheels Gravitational potential energy of the left front suspension, wheel and tire Gravitational potential energy for the right front suspension, wheel and tire Generalized forces associated with the rear suspension dampers, springs Generalized forces associated with the left front suspension damper and spring Generalized forces associated with the right front suspension damper and spring Generalized forces associated with all tire forces Generalized forces associated with the normalization constraint on the Euler parameters. Generalized constraint forces associated with the rear panhard rod and the trailing links. Generalized constraint forces associated with steering linkage attachments to the front suspension Generalized constraint forces associated with the upper and lower control arms for the front suspension

The first group of terms gives rise to scalar expressions in terms of qk , q&k and q&&k . The remaining terms give rise to scalar expressions containing only qk and q&k . The equations of motion are coupled in the q&&k terms but it is possible to solve for the generalized accelerations algebraically. Note that to simulate the motion of the vehicle using numerical

45

Table 1.4 - Organization of the Vehicle Model Derivations

Chapter

Vehicle Component

EOM Terms

2 3

Sprung Mass Front Suspension and Wheels

TSM, VSM, QC,k

4 5 6 7 8

Three Link Rear Suspension and Wheels Steering System Road Model Tire Model Driver Model

TLF, TRF, TFT, VLF, VRF, QLF,k, QRF,k, QCA,k TRS, TRT,VRS QTR,k QT,k

integration, the system of equations must be formulated and solved each time the integrator requests a function evaluation. The detailed derivations of the terms in the preceding equation are provided in the following chapters. As mentioned before, the terms are grouped by function. Table 1.4 shows the organization of the derivations.

46

2 Equations of Motion - Sprung Mass 2.1 Introduction The motion of the vehicle’s sprung mass is expressed in terms of seven generalized coordinates: three Cartesian displacements and four Euler parameters. A normalization constraint on the Euler parameters produces the desired six degree of freedom system. The terms in the equations of motion associated with the sprung mass are found by expressing the position of the sprung mass in terms of these generalized coordinates, differentiating to obtain the velocity of the sprung mass, finding the angular velocity of the mass in terms of the Euler parameters, formulating the kinetic energy and the potential energy of the body, and finally, differentiating the energy expressions to obtain the relevant terms in the equations of motion. A partial schematic of the vehicle is shown in Figure 2.1. The coordinate system E represents an Earth fixed coordinate system which is assumed to be inertial. The SM coordinate system is rigidly attached to the vehicle sprung mass and its origin is at the center of mass of the sprung mass. The unit vectors of the SM coordinate system are oriented with the s$1 axis pointing forward, the s$2 axis pointing out of the driver’s side of the vehicle and the

s$3 axis

pointing upwards. The s$1 − s$3 plane is chosen so that it is parallel to the symmetry plane of

47

s$3

s$2

SM e$3

s$1 rSM/E e$2

E

e$1 Figure 2.1: Earth Fixed and Vehicle Sprung Mass Coordinate Systems

the sprung mass; this is done to allow simplification of the steering model which is discussed in a later chapter. The position of the origin of the SM coordinate system with respect to the origin of the inertial coordinate system E is given by the vector rSM/ E .

2.2 Sprung Mass Kinetic and Potential Energy Terms Since the position of the vehicle is typically expressed with respect to the roadway surface and the data describing the roadway surface is expressed in terms of the Earth fixed coordinate system, it is desirable to represent the position of the vehicle in terms unit vectors associated with the E coordinate system: E

rSM/E = x e$1 + y e$2 + z e$3

(2.2.1)

The superscript at the upper left of the r indicates that the vector rSM/ E is written in terms of the unit vectors of coordinate system E. The subscript is intended to be read as “SM with respect to E” so E rSM/E is the position of point SM with respect to point E written in

48

terms of the unit vectors of the E coordinate system. The three degrees of freedom x, y and z represent the position of the sprung mass. Since the E system is inertial the velocity of the sprung mass can be obtained by direct differentiation giving E

vSM /E = E r&SM/ E = x& e$1 + y& e$2 + z& e$3

(2.2.2)

The angular orientation of the sprung mass can be described in several ways. In the past models of this type have ignored some of the rotational degrees of freedom which leads to an essentially one dimensional or two dimensional description for the angular position and angular velocity. This approach typically leads to a simple form for the equations of motion. While this approach can yield closed form solutions it compromises the accuracy of the model. In order to capture the subtleties of the motion it is necessary to utilize a representation which encompasses all three rotational degrees of freedom. This can be achieved by representing the orientation of the SM coordinate system with respect to the E system in terms of a rotation matrix. Several forms of rotation matrices have been used in the past with the most common ones being those based on Euler angle sequences. While this approach works well there are problems with the rotation matrix becoming singular at certain orientations of the body which can make solution of the equations of motion difficult. A better solution is to express the rotation matrix in terms of Euler Parameters. 1 The rotation matrix is shown below.

1

See Chapter 6 of Nikravesh’s text for a detailed explanation of Euler Parameters.

49

[E CSM ] = [G ][L]T 2 2 1  βSM βSM,1βSM, 2 − βSM, 0 βSM,3 βSM,1βSM, 3 + βSM, 0 βSM, 2  ,0 + βSM,1 − 2   (2.2.3) 2 2 1 = 2 βSM,1βSM, 2 + βSM, 0 βSM, 3 βSM, βSM, 2 βSM,3 − βSM, 0 βSM,1  0 + βSM, 2 − 2 2 2 1 β β  βSM, 0 + βSM,3 − 2  SM,1 SM,3 − βSM, 0 βSM, 2 βSM, 2 βSM,3 + βSM, 0βSM,1 

− βSM,1  [G ] = − βSM, 2  − βSM,3 − βSM,1  [L] = − βSM, 2  − βSM,3

βSM, 0

− βSM,3

βSM,3

βSM, 0

− βSM, 2

βSM,1

βSM, 0

βSM,3

− βSM,3

βSM, 0

βSM, 2

− βSM,1

βSM, 2   − βSM,1  βSM, 0   − βSM, 2   βSM,1  βSM, 0  

(2.2.4)

Note that there are four Euler parameters (βSM,i). In order for the rotation matrix to be normal a constraint on the Euler parameters is required. 2 2 2 2 βSM, 0 + βSM,1 + βSM, 2 + βSM, 3 = 1

The normalization constraint reduces the number of rotational degrees of freedom from four to the expected three degrees of freedom. The angular velocity of the sprung mass is written in terms of the Euler parameters as follows.

SM

ω SM/E

− β1  = 2[L]β& = 2 − β2  − β3

β0 − β3 β2

β3 β0 − β1

β&  − β2   &0   β  β1   &1  β2 β0   &   β3 

(2.2.1)

It is important that the angular velocity be written in terms of the sprung mass centroidal coordinate system (SM) so that the inertia tensor is constant with respect to the frame of reference.

50

With these considerations in mind the kinetic energy can be written as TSM =

1 1 T msm ( SM v TSM/E )( SM v SM/E ) + SM ω SM/E [SM J SM ] SM ω SM/ E 2 2

(2.2.2)

The velocity of the sprung mass written in terms of the body fixed coordinate system is related to the velocity in the inertial coordinate system by a simple coordinate transformation matrix: SM

vSM /E = [ E CSM ]T E v SM/E

β20 + β12 − 12  [E CSM ] = 2 β1β2 + β0β3  β1β3 − β0β2

β1β2 − β0β3 β1β3 + β0β2   β20 + β22 − 12 β2β3 − β0β1  β2β3 + β0β1 β20 + β23 − 12  

Substituting this result and the definition of

E

(2.2.3)

vSM/E into the expression for the kinetic

energy and simplifying gives

1 1 T mSM ( E v TSM/E [ E CSM ])([E CSM ]T E v SM/E ) + SM ω SM/E [SM J SM ]SM ω SM/E 2 2 1 1 T T = mSM ( E v SM/E )( E v SM/E ) + SM ω SM/E [SM J SM ]SM ω SM/E 2 2

TSM =

(2.2.4)

Given the form of the kinetic energy equation it is possible to find the terms in the equations of motion which are derived from the kinetic energy of the sprung mass. The relevant terms in the equations of motion are given by the expression

ETSM ,qk =

d  ∂TSM  ∂TSM  − dt  ∂q&k  ∂qk

(2.2.5)

Rather than expanding the expression for TSM and differentiating directly it is preferable to differentiate TSM and then substitute the derivatives of the linear velocity and the angular velocity into the result.

51

T

 ∂E vSM /E  ∂TSM 1 = mSM   2 ∂qk  ∂qk 

(

E

E 1 T  ∂ v SM /E  E v SM/E ) + mSM ( v SM/E )   2  ∂qk 

T

 ∂SM ω SM /E  1 ∂SM ω SM/ E  SM 1 SM T SM SM +   [ J SM ]( ω SM/ E ) + ( ω SM/E ) [ J SM ]  2  ∂qk 2   ∂qk 

(2.2.6)

The resulting terms are scalars which allows them to be transposed without affecting the equation. Note that the inertia tensor is symmetric so it is also unaffected by transposition. Thus, the first and second terms are equal and that the third and fourth terms are equal. E ∂TSM T  ∂ v SM/E  = mSM ( E vSM/ E )  + ∂qk  ∂qk 

(

SM

∂SM ω SM /E  T ω SM/ E ) [SM J SM ]   ∂qk 

(2.2.7)

Differentiation with respect to q&k leads to a result identical in form. To complete the Lagrangian the preceding result needs to be differentiated with respect to time (with q&k substituted for qk).

d ∂TSM  d    = mSM dt  ∂q&k  dt 

(

E

E T ∂ v  SM/ E v SM/ E  +  ∂q&k 

)

T

(

SM

 ∂E vSM / E   d E  v SM/ E  + mSM = mSM      ∂q&k   dt

T ∂SM ω SM/ E  ω SM / E [SM J SM ]   ∂q&k 

)

(

T

E

v SM/ E

)

T

d ∂E v SM/ E    dt  ∂q&k 

(2.2.8)

 ∂ ω SM/ E  SM   d SM ω SM / E  +   [ J SM ]   dt  ∂q&k  SM

+

(

SM

ω SM/ E

)

d  ∂SM ω SM / E  [ J SM ]   dt  ∂q&k 

T SM

Substituting the generalized coordinates x, y and z and their derivatives for qk and q&k into the preceding equations gives the following results:

52

ETSM ,x =

d ∂TSM  ∂TSM  − ∂x dt  ∂x& 

T E E T d  ∂ v SM /E  d E   ∂ vSM/ E  E v = mSM     + mSM ( vSM/E )    dt SM/ E   ∂x&  dt  ∂x& 

(2.2.9)

T  ∂ v SM/E  − mSM ( E vSM/ E )    ∂x  = mSM &x& E

ETSM , y =

d  ∂TSM  ∂TSM = mSM &y&  − dt  ∂y&  ∂y

(2.2.10)

ETSM ,z =

d ∂TSM  ∂TSM = mSM &z&  − dt  ∂z&  ∂z

(2.2.11)

Considering the rotational coordinates and applying Equations A.10 and A.11 gives E

ETSM ,βi =

d  ∂TSM  ∂TSM  − dt  ∂β&i  ∂βi T

 ∂SM ω SM / E  SM  E d SM  ω SM / E  =  [ J SM ] &  dt   ∂βi  +

(

SM

The derivatives of

(2.2.12)

 E d ∂SM ω SM/ E   ∂SM ω SM / E  T ω SM/ E [SM J SM ] −    dt    ∂β&i   ∂βi  

SM

)

ω SM / E are calculated in Appendix A. The first term in the preceding

expression is the only one containing second derivatives of the Euler parameters (due to the

d SM dt

&& . The remaining terms are nonlinear functions ω SM / E term) and it is linear in the β i

of the Euler parameters and their first derivatives. The potential energy of the sprung mass consists only of a gravitational potential energy term. The terms associated with the springs and dampers which support the sprung mass are included in the appropriate generalized force terms. This is done so that nonlinear

53

springs and dampers can be implemented without significantly revising the model. The mass terms associated with the suspension and wheels are lumped into the suspension potential energy terms and will be derived in a later section. It is assumed that gravity acts parallel to the e$3 axis. Given these considerations the potential energy for the sprung mass is simply

VSM = mSM gz

(2.2.13)

The terms in the equations of motion associated with the sprung mass potential energy are

EVSM ,qk =

∂VSM ∂qk

(2.2.14)

Clearly the only nonzero derivative is the one associated with the z degree of freedom. EVSM ,z = mSM g

(2.2.15)

2.3 Euler Parameter Constraints The normalization constraint on the Euler parameters used to represent the angular orientation must be incorporated into the equations of motion in some manner. While it is possible to eliminate one of the βi from the equations of motion using the constraint equation this leads to complicated nonlinear equations and the symmetry of the transformation matrix is destroyed. Although it is not as efficient, it is more elegant, to determine a set of constraint forces which can be applied to the model which will ensure that the constraint equation is satisfied at all times. The method of Lagrange multipliers

54

will be used to determine the constraint forces 2. The constraint forces are given by the equation QC ,βi = λC a βi

(2.3.1)

where

a βi =

∂ 2 β0 + β12 + β22 + β23 − 1) ( ∂βi

(2.3.2)

Evaluating Equation 3.39 gives a β0 = 2β0

a β1 = 2β1

a β2 = 2β2

a β3 = 2β3

(2.3.3)

Note that there is now one additional unknown λC (the constraint force) which must be determined. It will be necessary to determine the value of λC for each integration function evaluation. The λC s can be determined by appending them to the acceleration vector and augmenting the equations of motion with the second derivative of the constraint equations. The first and second derivatives of the normalization constraint equation are as follows.

d 2 β0 + β12 + β22 + β32 − 1 = 0} { dt ⇒ β0β&0 + β1β&1 + β2β&2 + β3β&3 = 0

{

}

d β β& + β1β&1 + β2β&2 + β3β&3 = 0 dt 0 0 && + β β && && && &2 &2 &2 &2 ⇒ β0 β 0 1 1 + β2 β2 + β3β3 + β0 + β1 + β2 + β3 = 0

2

See Meirovitch’s “Methods of Analytical Dynamics”, Chapter 2 for derivation of the equations.

55

(2.3.4)

(2.3.5)

3 Equations of Motion - Front Suspension and Wheels 3.1 Introduction The relatively complicated geometry of the front suspension creates a dilemma. To obtain an accurate model it is necessary to represent the linkages between the spindle and the sprung mass as exactly as possible. The best way of doing this is to allow the spindle to have six degrees of freedom and then limit its motion via constraint equations. This type of model is not optimal if one is concerned about computational efficiency. The opposing approach is to assign two degrees of freedom (steering angle and spring length for instance) to the spindle. This approach leads to a potentially more efficient model but the difficulty in deriving the appropriate equations of motion is prohibitive. For this reason the first modeling technique discussed is chosen in spite of computational costs. A schematic of the sprung mass and the upper and lower control arms on the left side of the vehicle is shown in Figure 3.1. The spindle is omitted for clarity. The upper and lower control arms are located with respect to the sprung mass by the vectors SM

rUC/SM and

SM

rLC/SM respectively. Each control arm has a coordinate system associated

with it. The UC coordinate system is associated with the upper control arm and has the unit vectors u$ 1 = SM u$ UA and u$ 2 . The vector

SM

u$ UA is a constant unit vector (with respect

to the sprung mass coordinate system) which points along the direction of the rotation axis

56

s$ 3 s$1 u$ 1 SM

UC

u$ UA

SM

$l LA

$l 1

u$ 2

SM

rUC/SM

LC

SM

SM

rLC/SM

s$ 2

$l 2 Figure 3.1: Schematic Showing the Front of the Sprung Mass and the Control Arms

of the control arm. The origin of the coordinate system is located such that the u$ 2 unit vector points at the ball joint where the control arm connects to the spindle. The coordinate system for the lower control arm (LC) is set up analogously with unit vectors $l = SM $l and $l . The control arms are assumed to be massless. This avoids the need for 1 UA 2

the control arms to be modeled as separate bodies which would increase the number of degrees of freedom of the model. The mass of the control arms can be accounted for by distributing it to the spindle and the sprung mass. A schematic showing the spindle and the control arms is shown in Figure 3.2. The spindle has a coordinate system attached at it’s mass center (SP) with unit vectors p$ i . The upper and lower ball joints (UJ and LJ respectively) which connect the spindle to the control arms are located by the vectors

SP

57

rUJ /SP and

SP

rLJ /SP . The coordinate system

associated with the right spindle is

UC p$ 3

oriented in the same fashion as the

u$ 1

coordinate system for the left spindle (i.e. both y-axes point to the driver’s

SP

SP

p$ 1

u$ 2

rUJ /SP

left and both x-axes point straight SP

ahead). The p$ 2 axis is parallel to the

p$ 2

axis of rotation of the wheel and tire. This is done so that the inertia tensor

rLJ /SP

LC

$l 1

$l 2

associated with the wheel, tire and Figure 3.2: Schematic of the Spindle and the Control Arms

brake disc (or drum) retains its symmetry and its constant value with respect to the SP coordinate system.

The spring and damper are assumed to be connected to the lower control arm. The location of the point of attachment is specified relative to the LC coordinate system by the vector

LC

rLM / LC . The subscript LM indicates the lower mount. The upper mount for the

spring and damper is assumed to be located on the sprung mass at a location given by the constant vector

SM

rUM /SM . The subscript UM indicates the upper mount. The forces

produced by the front springs and dampers are applied to the model via generalized forces. The front suspension is modeled using a 7 degree of freedom spindle whose motion is constrained. There are a total of five constraint equations: four physical constraints which represent the control arms and one normalization constraint on the Euler

58

parameters. This leaves two unconstrained degrees of freedom which can be roughly equated with steer angle and vertical position of the wheel with respect to the sprung mass. The steer angle degree of freedom is eliminated by another constraint equation associated with the steering system model which is discussed in a later chapter.

3.2 Front Spindle Kinetic and Potential Energy Terms The terms in the differential equations associated with the motion of the spindle are calculated below. The kinetic energy terms associated with the linear motion of the spindle include the mass of the wheel and tire assembly. The rotational kinetic energy of the spindle, but not of the wheel and tire assembly, is also included. The position of the spindle relative to the inertial coordinate system is given by the vector E

r SP/ E = xSP e$ 1 + ySP e$ 2 + zSP e$ 3

(3.2.1)

The angular orientation is given by the four Euler parameters βSP,i . The form of the kinetic and potential energy expression for the spindle are identical to the formulations for the sprung mass. The resulting terms in the equations of motion are identical in form as well. ETSP ,x = mSP &x&SP

(3.2.2)

ETSP , y = mSP &y&SP

(3.2.3)

ETSP ,z = mSP &z&SP

(3.2.4)

59

T

ETSP ,βSP,i

∂SP ω SP/ E  SP  d SP   ω SP/ E  =  ∂β&  [ J SP ]  dt    SP, i +

(

SP

ω SP/ E

 E d  ∂SP ω   ∂SP ω SP/ E  T SP SP/ E   [ J SP ]  −   dt  & ∂ β ∂ β     SP,i   SP,i

(3.2.5)

)

The terms in the equations of motion associated with the spindle potential energy are EVSP ,z = mSP g

(3.2.6)

3.3 Front Wheel and Tire Rotational Energy Terms The front wheels are assumed to be rigidly affixed to the spindle/steering knuckle assembly and to rotate about an axis parallel to the p$ 2 axis of the spindle coordinate system. The terms in the equations of motion associated with the kinetic energy of the front wheels due to rotational motion are derived here. Due to the rotational symmetry of the wheel assembly and the alignment of the p$ 2 axis with the axis of rotation it is not necessary to utilize a distinct coordinate system for the wheels; the inertia tensor is constant in the SP coordinate system. The kinetic energy for one of the wheels can be written TWH =

1 SP T ω wheel / E [SP J wheel ]SP ω wheel / E 2

(3.3.1)

The angular velocity vector of the wheel with respect to the SP coordinate system has constant direction but variable magnitude. The total angular velocity of the wheel includes the angular velocity of the SP coordinate system. SP

& ω wheel / E = SP ω SP/ E + φ wheel

60

SP

p$ 2

(3.3.2)

where φwheel is the angular velocity degree of freedom (scalar) associated with the rotation of the wheel. The value is time dependent and is dictated by the interaction of the vehicle model and the tire model. The terms in the equations of motion related to the rotational motion of the wheel are given by E

ETWH ,qk =

d  ∂TWH  ∂TWH  − dt  ∂q&k  ∂qk

(3.3.3)

Substituting for the kinetic energy and differentiating gives T

ETWH ,qk

 E d SP   ∂SP ω wheel / E  SP = [ J ] ω wheel / E   wheel   ∂q&k   dt  +

(

SP

ω wheel / E

)

 E d  ∂SP ω wheel / E  ∂SP ω wheel / E  [ J wheel ]   − ∂qk   dt  ∂q&k 

(3.3.4)

T SP

The only degrees of freedom which generate non-zero results for the expression above are the βSP,i and φwheel. The derivatives of the wheel angular velocities are calculated as follows ∂SP ω wheel / E ∂SP ω SP/ E = ∂βSP,i ∂βSP,i

∂SP ω wheel / E ∂SP ω SP/ E = ∂β&SP,i ∂β&SP,i

∂SP ω wheel / E =0 ∂φwheel E

E d SP d ω wheel / E = dt dt

E

(

SP

(3.3.5)

∂SP ω wheel / E SP = p$ 2 & ∂φ

(3.3.6)

wheel

)

&& ω SP/ E + φ wheel

SP

p$ 2 +

SP

& ω SP/ E ×φ wheel

d  ∂SP ω wheel / E  E d ∂SP ω SP/ E      =   dt   ∂β&SP,i  dt  ∂β&SP,i 

61

SP

p$ 2

(3.3.7)

(3.3.8)

E

d ∂SP ω wheel / E d SP d = p$ 2 = & dt ∂φwheel dt dt E

SP

(

SP

)

p$ 2 +

SP

ω SP/ E ×SP p$ 2

(3.3.9)

= SP ω SP/ E ×SP p$ 2

For qk = βSP,i and substituting for the derivatives of

SP

ω wheel / E , using the expressions

above, the result becomes  E d  ∂SP ω T  ∂SP ω T  SP SP/ E SP/ E    ETWH ,βSP,i =  − [ J wheel ]   dt  & ∂βSP,i   ∂βSP,i   

(

T  ∂SP ω SP/  SP E d E  J [ ] +  wheel   ∂β&   dt  SP,i 

(

SP

ω SP/ E

)

SP

SP & p$ 2 ω SP/ E + φ wheel

)

 SP && SP p$ + SP ω & p$ 2  + φ 2 wheel SP/ E × φwheel 

(3.3.10)

For qk = φwheel the expression simplifies as follows

ETWH ,φ =

(

+

(

SP

(

)

T

ω SP/ E ×SP p$ 2 [SP J wheel ]

SP

 SP d p$ T2 [SP J wheel ]  dt

)

(

SP

SP

SP & ω SP/ E + φ p$ 2 wheel

)

& SP p$ + ω SP/ E + &φ wheel 2

)

SP

 SP & ω SP/ E × φ p$ 2  wheel 

(3.3.11)

3.4 Generalized Forces for Springs and Dampers To calculate the generalized forces associated wi th the front suspension it is necessary to calculate the virtual work done by the spring and the damper. The virtual work done by the spring or damper is equal to the force exerted by the spring or damper multiplied by the change in length of the spring or damper due to a virtual displacement of the system: δW = E FSD ⋅δE LSD E

(3.4.1)

FSD is the force vector due to the spring or damper. The vector δE LSD is the change in

length due to a virtual displacement of the upper and lower spring mounts. To determine

62

the length vector is it necessary to find the position of the upper and lower mounts in terms of the generalized coordinates. Since the spring and damper are typically attached to the lower control arm for type of vehicle being considered the derivation below is carried out using the lower control arm vectors defined in Figure 3.1. The unit vector

E

$l is given in terms of the 1

sprung mass coordinate system as part of the vehicle model specification and is determined trivially. The

E

$l unit vector is found by determining the vector from the origin of the 2

control arm coordinate system to the ball joint on the spindle ( E l BJ / CA ), subtracting those E

components of the resulting vector which are parallel E

E

l BJ / CA =

(

E

r SP/ E +

E

$l = 2

$l = [ C ]SM $l 1 E SM LA

[C ]

SP

E

SP

(3.4.2)

r LJ /SP − E r SM / E −

( −(

) $l $l ) $l

l BJ / CA −

E T E E BJ /CA 1 1

E

l BJ / CA

E T E E BJ /CA 1 1

E

l

l BJ / CA −

l

(

SM

E

SM

r LC/SM

)

$l

) $l

(3.4.3)

(3.4.4)

E T E E BJ /CA 1 1

l

−

E T E BJ / CA BJ /CA

l

[C ]

$l

E

=

$l and normalizing the result. 1

l

(

$l

E T E BJ / CA 1

l

)

2

Using these definitions the position of the lower spring mount can be written as E

r LM / E = E r SM / E +

where E rLM / LC = LC x LM / LC E $l1 +

LC

[C ]

SM

E

SM

y LM / LC E $l 2 +

LC

r LC/SM + E rLM / LC

(3.4.5)

z LM / LC E $l 3 is a constant vector (in the control

arm coordinate system) which locates the spring or damper mount with respect to the origin of the LC coordinate system. To simplify the development of the equations it is

63

assumed that the lower mounts are in the plane of the control arm (i.e. Note that

SM

LC

z LM / LC = 0 ).

r LC/SM is a constant vector with respect to the SM coordinate system. The

upper mount is typically attached directly to the sprung mass so its position can be written as E

r UM / E = E r SM / E +

[C ]

SM

E

SM

r UM/SM

(3.4.6)

The length of the spring or the damper is then E

LSD = E r UM / E − E r LM / E

[C ] ]( r

= E r SM / E + = [E CSM

SM

E

SM

SM

UM /SM

r UM /SM − E r SM / E −

−

SM

r LC/SM )−

LC

[C ]

SM

E

SM

x LM / LC E $l1 −

r LC/SM − E rLM / LC

LC

(3.4.7)

y LM / LC E $l 2

Differentiating to obtain the virtual displacement gives

 ∂[E CSM ] SM  δE LSD = ∑   ( r UM /SM − i  ∂βSM, i  −

LC

(

)

x LM / LC δE $l1 −

LC

(

SM

r LC/SM )δβSM, i

y LM / LC δE $l 2

(3.4.8)

)

The derivatives of the unit vectors of the lower control arm coordinate system are

δE $l1 =

∑ i

∂[E CSM ]SM   $l LA δβSM, i  ∂β   SM, i 

64

(3.4.9)

$l

) (

E

$l

)

(

δ $l 2 = δE l TBJ / CA E l BJ / CA −  

(

E T E BJ /CA 1

δE $l 2 = − E l TBJ / CA E l BJ / CA −  

(

E T E BJ / CA 1

E

(( l

⋅

E T BJ /CA

(

−

l

$l

E T E BJ / CA 1

l

(

+ E l TBJ / CA E l BJ / CA −  

(

l

2

) $l )δl E

−

  

2

−

1 2

3

$l

E T E 1 BJ / CA

l

−

BJ / CA

)

2

(

−

(

BJ /CA

1

(

(

⋅ E l TBJ / CA E $l1

)(

$l

E T E 1 BJ /CA

l

E

  

2

E

BJ / CA

 −  E l$2TδE l BJ /CA 

(

1

$l

l

)

(

1 2

$l

)l

δE $l1

l

δE $l1

)

) )

$l δE $l 1

l

E

BJ / CA

E T E BJ / CA BJ / CA

E T BJ / CA

) $l )

l

E T E 1 BJ / CA

E T E BJ / CA 1

1

) $l )

E T E E BJ / CA 1 1

) {δl $l )δ $l δ $l ) $l − ( l )−  l l − ( l $l ) 

δE $l 2 = E l TBJ / CA E l BJ / CA −   E$ l + ElT − δE l T

(

l

1 2

   −

$l

E T E E BJ / CA 1 1

E T E E T 1 BJ / CA BJ /CA

⋅ δE l BJ / CA − δE l TBJ / CA E $l1 + E l TBJ /CA δE $l1 E $l1 −

(

(

 2 ⋅ El BJ / CA −  

E

1

l BJ / CA −

E

1

E T E BJ / CA 1

2

−

1 2

)) $l } E

(3.4.10)

2

Simplifying a bit further gives the final result.

{

 1  E E T E E T E E δE $l 2 =   δ l BJ /CA − δ l BJ / CA $l1 + l BJ / CA δ $l1 $l1 −  mag 

(

)

  1  E T E −  E $l 2TδE l BJ /CA −   l BJ / CA $l1  mag  

(

)

mag ≡ E l TBJ / CA E l BJ / CA −  

(

(

$l

E T E 1 BJ / CA

l

)

2

−

  

)(

E T BJ /CA

l

(

l

)

$l δE $l 1

E T E 1 BJ /CA

   δE $l1 E l$2    

)

(3.4.11)

1 2

where

δE l BJ /CA = (1 0 0)δxSP + ( 0 1 0)δ ySP + ( 0 0 1)δ z SP

− (1 0 0)δ xSM − ( 0 1 0)δ ySM − ( 0 0 1)δ z SM +

∑ i

∂[E CSP ]SP r LJ /SP δβSP,i − ∂βSP,i

∑ i

(3.4.12)

∂[E CSM ]SM r LC/SM δ βSM,i ∂βSM,i

Substitution of the expression for δE l BJ / CA into the expression for δE $l 2 produces an unwieldy result. Rather than write the complete expression it is preferable to break it 65

down into components based on the degree of freedom. The complete expression for δE $l 2 is obtained by summing the individual terms.

{

( $l (1

0 0) E $l1 −

{

( $l ( 0

1 0) E $l1 −

{

( $l ( 0

0 1) E $l1 −

{

( $l (1

0 0) E $l1 −

{

( $l ( 0

1 0) E $l1 −

{

( $l ( 0

0 1) E $l1 −

 1  δxSM ⇒ −   (1 0 0) −  mag   1  δySM ⇒ −   ( 0 1 0) −  mag   1  δzSM ⇒ −   ( 0 0 1) −  mag   1  δxSP ⇒   (1 0 0) −  mag 

E T 1

E T 1

E T 1

E T 1

 1  δySP ⇒   ( 0 1 0) −  mag 

E T 1

 1  δzSP ⇒   ( 0 0 1) −  mag 

E T 1

α ≡−

∂βSM,i

r LC/SM , γ≡

 1  δβSP,i ⇒   α−  mag 

{ (

E

)

) }

(3.4.15)

) }

(3.4.16)

) }

(3.4.17)

) }

(3.4.18)

1 0) E $l 2 δySM

)

( $l ( 0

)

( $l (1

0 0) E $l 2 δxSP

)

( $l ( 0

1 0) E $l 2 δySP

)

( $l ( 0

0 1) E $l 2 δzSP

E T 2

E T 2

E T 2

E T 2

)(

E T 2

(

l

E T BJ / CA

l

∂βSM,i

(

)

$l γ

E T E 1 BJ / CA

∂[E CSM ]SM

$l Tα E $l − 1 1

∂[E CSP ]SP α≡ r LJ /SP ∂βSP,i

0 1) E $l 2 δzSM

( $l ( 0

(

∂[E CSM ]SM

(3.4.14)

)

E T 2

)

)

) }

0 0) E $l 2 δxSM

  1 E T −  E l$2Tα −   l BJ / CA E $l1 mag   

(

(3.4.13)

( $l (1

 1  δβSM,i ⇒   α − α T E $l1 + E l TBJ / CA γE $l1 −  mag 

{ (

) }

)

  γ)E l$2 δ βSM,i  

(3.4.19)

$l LA

) }

$l α E $l δβ 2 SP,i

E T 2

(3.4.20)

66

The final step in determining the virtual work is calculating the magnitude and direction of the force vector. The force exerted by the spring is typically a function (possibly nonlinear) of the length. The force exerted by the damper is a potentially nonlinear function of the time derivative of the length. The length vector for the spring or damper was already determined to be E

LSD = [E CSM ]( SM r UM /SM −

SM

r LC/SM )−

LC

x LM / LC E $l1 −

LC

y LM / LC E $l 2

(3.4.21)

& x LM / LC E $l1 −

LC

& y LM / LC E $l 2

(3.4.22)

Taking the time derivative of this expression gives E

SM & =[ C & L r UM/SM − SD E SM ](

SM

r LC/SM )−

LC

where E

E

{

1  E& & $l =  −  l  2  mag  BJ / CA

&$ & ]SM $l l1 = [E C SM LA

( &l

)

E& $l + E l T $ E $l − 1 BJ / CA l1

E T E 1 BJ /CA

  1  E T −  E $l 2T E &l BJ / CA −   l BJ/ CA E $l1  mag  

(

mag ≡ E

Once

E

)

[l

E T E BJ / CA BJ / CA

l

&l BJ / CA =

LSD and

E

(

E

r&SP/ E +

(

(

−

)]

2 − E T E 1 BJ / CA

$l

l

[ C& ]

SP

E

(3.4.23)

SP

)(

(

$l

) &$l

E T E E 1 1 BJ / CA

l

)

 &$ E $  l  l2    

E T E 1 BJ / CA

l

(3.4.24)

1 2

r LJ /SP − E r&SM / E −

[ C& ]

SM

E

SM

r LC/SM

)

(3.4.25)

& are known the force generated by the spring or damper may be L SD

calculated. In most cases the length of the spring is not given directly by

E

LSD . An

example of this is a coil-over shock assembly. For this type of assembly the preload on the spring is set by adjusting a ring which moves up and down the body of the damper on a

67

thread. For modeling purposes the preload on the spring can be specified by providing the distance along the coil over shock assembly from the appropriate mount to the end of the spring as an input. If the force-displacement curve of the spring is modeled as a polynomial containing linear and cubic terms the expression for the spring force can be written as

∆ L = fS − E

FS

(

E

(

E

LS − pS

)(

)

)

$ LS = k 0 ∆ L + k1∆ L3 E L S

(3.4.26)

where fs is the free length of the spring, ps is the preload length as described above, ELS is the magnitude of the vector ELS,

E

$ is the unit vector along the direction of ELS, ∆L is L S

the amount of compression applied to the spring, and k0 and k1 are the linear and cubic stiffness coefficients. The force generated by the damper is typically dependent on the rate of change of the length of the damper. The amount of damping is usually dependent on whether the damper is extending (rebound) or compressing (jounce). The force generated by the damper can be modeled as follows: E

FD

( L& )= − (c E

D

E R0

)

3 $ L&D + cR1 E L&D E L D

(3.4.27)

where E

d L&D = dt

[L E

T E D

]= 1

LD

2

E

$T EL & L D D

$ is the unit vector along the direction of ELD. and E L D

68

(3.4.28)

The virtual work δW = E FSD ⋅δE LSD is obtained by combining the above results. The generalized forces are the coefficients of the corresponding virtual displacements. The results are summarized in Table 3.3 at the end of the chapter.

3.5 Constraint Forces for the Control Arms The motion of the front spindle is constrained by the action of the control arms and the action of the steering linkage. The control arms and steering linkage are modeled using constraint equations. The discussion of the steering system constraints is covered in the chapter on the steering model. The constraint equations for the control arms are derived below. There are a total of four control arms (two per side) used in the typical SLA front suspension used on most race cars. The derivation procedure is identical for each of the control arms. To eliminate redundancy the following derivation is done for a generic control arm. The vectors describing the generic control arm are depicted in Figure 3.3. Each control arm can be represented by two constraint equations. The first constraint forces the ball joint attachment point (BJ) on the spindle to remain a fixed distance (i.e. the length of the control arm LCA) from the origin of the control arm coordinate system (CA) which is affixed to the sprung mass. This can be written mathematically as

69

E

r BJ / CA = E r SP/ E +

[C ]

SP

E

SP

[r E

T E BJ / CA

r BJ /SP − E r SM / E −

[C ]

SM

E

SM

r CA /SM

(3.5.1)

r BJ /CA ] − LCA = 0 1 2

(3.5.2)

The second condition is that the same vector be perpendicular to the axis of rotation of the control arm. This can be expressed mathematically by noting that the inner product of the vectors must be zero: E

rBJT /CA [E CSM ]SM c$ CAA = 0

(3.5.3)

The generalized constraint forces are of the form QCA,qi = λCA a qi

(3.5.4)

where

a qi =

∂E T E  rBJ /CA r BJ / CA ∂qi 

[

]− 1 2

 ∂  = E r$BJT /CA  E r BJ /CA  LCA   ∂q 

(3.5.5)

i

or

s$ 3 s$1

p$ 3

SM c$ 1 = c$ CAA

SM

CA p$ 1

SP SP

rBJ /SP

rCA / SM

s$ 2

BJ c$ 2

p$ 2 Figure 3.3: Schematic of a Generic Control Arm

70

a qi =

∂ ∂qi

(

E

rBJT /CA [E CSM ]SM c$ CAA

) (3.5.6)

 ∂   ∂  SM =  E rBJT /CA [E CSM ]SM c$ CAA + E rBJT / CA   c$ CAA [ E C SM ] ∂qi   ∂qi 

{

For qi ∈ xsp , ysp , zsp , xsm , ysm , z sm

∂ E r = ∂xsp BJ /CA

(1

}the

∂ E ∂qi BJ / CA

r

∂ E r = ∂ysp BJ /CA

0 0)

(0

term evaluates to

1 0)

∂ E r = ∂z sp BJ / CA

(0

0 1)

∂ E ∂ E ∂ E r BJ / CA = − (1 0 0) r BJ / CA = − ( 0 1 0) r = − ( 0 0 1) ∂xsm ∂y sm ∂z sm BJ / CA

{ }

∂ E ∂qi BJ / CA

For qi ∈ βSP,i the

r

term evaluates to

 SP  ∂ ∂ E  r BJ /SP r BJ / CA =  [ E C SP ] ∂βSP,i   ∂βSP,i

{

}

For qi ∈ βSM, i the

∂ E ∂qi BJ / CA

r

(3.5.7)

(3.5.8)

term evaluates to

 SM  ∂ ∂ E  r CA /SM r BJ /CA = −  [ E CSM ] ∂βSM, i   ∂βSM,i

(3.5.9)

In order to determine the Lagrange multipliers it is necessary to append the second derivative of the constraint equations to the equations of motion. The first and second derivatives are calculated below. For the length constraint on the control arm:

d dt ⇒

{[ r E

T E BJ / CA

[r E

T E BJ / CA

}

r BJ /CA ] − LCA = 0 1 2

r BJ/ CA ] −

1 2

E

rBJT / CA E r&BJ / CA = 0

(3.5.10)

⇒ E rBJT / CA E r&BJ /CA = 0 d E T E { rBJ / CA r&BJ / CA = 0} dt ⇒ E rBJT / CA E &r&BJ / CA + E r&BJT / CA E r&BJ /CA = 0

71

(3.5.11)

For the orthogonality constraint:

{

}

d E T rBJ / CA [E CSM ]SM c$ CAA = 0 dt & ]SM c$ = 0 ⇒ E r&BJT / CA [E CSM ]SM c$ CAA + E rBJT /CA [E C SM CAA

{

}

[ ] [ ]

d E T SM & r&BJ / CA [E CSM ]SM c$ CAA + E rBJT /CA E C c$ CAA = 0 SM dt T SM E T SM & $ & ⇒ E &r&BJ c$ CAA + E rBJT / CA / CA [ E CSM ] c CAA + 2 rBJ / CA E CSM

[

&& E C SM

(3.5.12)

]

SM

(3.5.13)

c$ CAA = 0

The derivatives of the vector E r BJ / LC are E

E

r&BJ / CA = E r&SP/ E +

[ C& ]

r BJ /SP − E r&SM / E −

[ C& ]

r CA /SM

(3.5.14)

r&&BJ / CA = E &r&SP/ E +

[ C&& ]

r BJ /SP − E &r&SM / E −

[ C&& ]

r CA /SM

(3.5.15)

SP

E

SP

SP

E

SP

72

SM

E

SM

SM

E

SM

3.6 Summary of Results Table 3.1 - Front Suspension Kinetic and Potential Energy Terms

Generalized Coordinate xSP ySP z SP

∂TSP q&k

d dt

∂TSP qk

+

∂VSP qk

)

Term

mSP x&&SP mSP &y&SP mSP &z&SP + mSP g

βSP,i

T

 ∂SP ω SP/ E  SP  d SP    ω SP/ E   ∂β&  [ J SP ]   dt   SP, i +

(

SP

 E d  ∂SP ω   ∂SP ω SP/ E  T SP/ E   ω SP/ E [SP J SP ]  −   dt   ∂β&SP,i   ∂βSP,i   

)

Table 3.2 - Wheel and Tire Rotational Energy Terms

Generalized Coordinate φSP

βSP,i

( ( )−

( ( )− ) d dt

∂TSP q&k

∂TSP qk

Term

(

SP

+

(

(

)

T

ω SP/E ×SP p$ 2 [SP J wheel ] SP

 SP d p$ T2 [SP J wheel ]  dt

)

(

SP

SP

SP & ω SP/E + φ p$ 2 wheel

)

& SP p$ + ω SP/ E + &φ wheel 2

 E d  ∂SP ω T  ∂SP ω T  SP SP/E SP/ E     [ J wheel ]  − ∂β  dt   & ∂ β  SP,i   SP,i 

(

T  ∂SP ω SP/E  SP E d   +  [ J wheel ]  dt  ∂β&SP,i 

(

73

SP

)

SP

) SP

 SP & ω SP/ E × φ p$ 2  wheel 

SP & p$ 2 ω SP/E + φ wheel

& SP p$ + ω SP/E + &φ 2 wheel

SP

)

 SP & p$ 2  ω SP/E × φ wheel 

(

Table 3.3 - Generalized Forces due to Spring or Damper QSD, k

Generalized Coordinate xSP , ySP , zSP

)

Generalized Force  LC y LM / LC E T −  F α − E $l1Tα E $l1 − E $l 2Tα E $l 2  mag  SD α = (1 0 0) for xSP , α = (0 1 0) for ySP and α = (0 0 1) for zSP

{ (

)

) }

(

 LC y LM / LC E T   F α − E $l1Tα E $l1 − E $l 2Tα E $l 2  mag  SD α = (1 0 0) for xSM , α = (0 1 0) for ySM and α = (0 0 1) for zSM

{ (

xSM , ySM , zSM

βSP,i

 LC y LM / LC E T −  F α−  mag  SD

{ (

α≡ βSM,i E

 ∂[E CSM ] T  FSD   ∂βSM,i 

(

SM

r UM /SM −

)

E

)

$l T α E $l − 1 1

SM

)

r LC/SM −

{ (

LC

)

  1  E T E −  E $l 2Tα −   l BJ / CA $l1  mag  

)

(

α ≡−

(

) }

$l α E $l 2

E T 2

∂[E CSP ]SP r LJ /SP ∂βSP,i  ∂[E CSM ]SM T  $l LA  x LM / LC E FSD ∂ β  SM, i 

 LC y LM / LC E T TE E T E −  FSD α − α $l1 + l BJ / CA γ $l1 −  mag 

(

) }

(

)(

E T BJ /CA

l

(

   γE $l 2    

l

)

$l γ

E T E 1 BJ / CA

)

∂[E CSM ]SM ∂[E CSM ]SM $l r LC/SM , γ≡ LA ∂βSM,i ∂βSM, i

74

(

Table 3.4 - Generalized Forces due to Control Arm Length Constraint QCA, k

Generalized Coordinate xSP

ySP z SP xSM ySM zSM

Generalized Force λCA E r$BJT / CA (1 0 0) λCA E r$BJT / CA ( 0 1 0) λCA E r$BJT / CA ( 0 0 1)

− λCA E r$BJT / CA (1 0 0) − λCA E r$BJT / CA (0 1 0) − λCA E r$BJT / CA (0 0 1)

βSP,i

 ∂[ C ] λCA E r$BJT / CA  E SP  SP r BJ /SP  ∂βSP,i 

βSM,i

 ∂[ C ] − λCA E r$BJT / CA  E SM  SM r CA /SM  ∂βSM,i 

75

)

(

Table 3.5 - Generalized Forces due to Control Arm Orthogonality Constraint QCA, k

Generalized Coordinate xSP

ySP z SP xSM

ySM

zSM βSP,i

βSM,i

)

Generalized Force λCA,2 (1 0 0)[E CSM ]SM c$ CAA λCA,2 ( 0 1 0)[E CSM ]SM c$ CAA λCA,2 ( 0 0 1)[E CSM ]SM c$ CAA

   ∂[E CSM ] SM  c$ CAA − (1 0 0)[E CSM ]SM c$ CAA  λCA,2 E rBJT / CA   ∂x      SM      ∂[E CSM ] SM  c$ CAA − ( 0 1 0)[E CSM ]SM c$ CAA  λCA,2 E rBJT / CA   ∂y     SM       ∂[E CSM ] SM  c$ CAA − ( 0 0 1)[E CSM ]SM c$ CAA  λCA,2 E rBJT / CA   ∂z      SM    ∂ E CSP  SP  r BJ /SP E CSM SM c$ CAA λCA,2   ∂β   SP, i 

[ ]

[

]

 ∂[E CSM ] SM  ∂[ C ]  r CA /SM [E CSM ]SM c$ CAA + E rBJT / CA  E SM  SM c$ CAA λCA,2   ∂β   ∂β    SM,i  SM,i 

76

4 Equations of Motion - Three Link Rear Suspension 4.1 Introduction The equations derived in this chapter are for a rear suspension geometry which consists of a solid rear axle supported by a three or four link mechanism (in the absence of bushing compliance the fourth link is redundant and needn’t be explicitly modeled). This geometry is similar to the layout used by the Legends series race cars and some dragsters. The motion of the rear suspension is expressed in terms of three position coordinates and four Euler parameters. As in the previous section the terms in the equations of motion associated with the rear unsprung mass are found by formulating the kinetic energy and the potential energy of the body and differentiating the energy expressions with respect to the generalized coordinates.

4.2 Unsprung Mass Kinetic and Potential Energy Terms In order to calculate the kinetic energy associated with the rear suspension of the vehicle model it is first necessary to calculate the velocity of the origin of the centroidal coordinate system associated with the rear suspension. The velocity is found by differentiating the position vector which is written as E

r RS/ E = x RS e$ 1 + y RS e$ 2 + z RS e$ 3

77

(4.2.1)

The velocity vector can be expressed in terms of several different coordinate systems. The Earth fixed inertial coordinate system (E) is chosen so that the expressions for velocity, position and acceleration have simple dependence on the variables x, y, z and their derivatives and to simplify computations of tire to road contact. Differentiating Equation 4.2.1 gives E

v RS/ E = ( x&RS

z&RS )

y&RS

(4.2.2)

The kinetic energy of the rear suspension is TRS =

1 1 mRS ( E v TRS/ E )( E v RS/ E ) + 2 2

RS

ω TRS/ E [RS J RS ]RS ω RS/ E

(4.2.3)

The mass of the rear suspension includes the mass of the rotating components (wheels, tires and brakes) as well as the mass of the rear axle, differential and suspension. The inertia tensor [RS J RS ] does not included the inertia of the rotating components. The terms related to the angular velocity of the wheels are considered separately below. The form of the kinetic energy for the rear suspension is identical to the expressions obtained for the sprung mass and the front unsprung masses. The resulting terms in the equations of motion are identical as well. ETRS ,x RS = mRS &x&RS

(4.2.4)

ETRS , yRS = mRS &y&RS

(4.2.5)

ETRS , z RS = mRS &z&RS

(4.2.6)

78

T

ETRS ,βRS,i

 ∂RS ω RS/ E  RS  d RS   ω RS/ E  =  ∂β&  [ J RS ]  dt    RS, i +

The derivatives of

(

RS

RS

ω RS/ E

 d ∂RS ω   ∂RS ω RS/ E  T RS RS/ E   [ J RS ]  −   dt  & ∂ β ∂ β    RS, i    RS,i

(4.2.7)

)

ω RS/ E which appear in the above expression are calculated in

Appendix A. The potential energy of the rear suspension consists only of a gravitational potential energy term. It is assumed that gravity acts parallel to the e$3 axis. Given that this is the case the potential energy is simply

VRS = mRS gz RS

(4.2.8)

The terms in the equation of motion associated with the rear suspension potential energy are

EVRS ,qk =

∂VRS ∂qk

(4.2.9)

Clearly the only nonzero derivative is the one associated with the zRS degree of freedom. EVRS ,z = mRS g

(4.2.10)

4.3 Rear Wheel Rotational Energy Terms The rear wheels are assumed to be rigidly affixed to the rear suspension and to rotate about the

RS

r$2 axis of the RS coordinate system. The kinetic energy and potential

energy terms associated with the translational motion of the rear wheels were derived by lumping mass of the wheels, tires and brakes with the rest of the rear suspension. The kinetic energy associated with the rotational motion of the wheels must be treated

79

separately due to the relative rotation between the wheels and the RS coordinate system. Due to the rotational symmetry of the wheels it is not necessary to utilize a distinct coordinate system for each wheel; the inertia tensor is constant in the RS coordinate system. The kinetic energy for a rear wheel can be written TRT =

1 2

RS

ω Twheel / E [RS J wheel ]RS ω wheel / E

(4.3.1)

The angular velocity vector associated with the wheel has constant direction but variable magnitude with respect to the RS coordinate system. The total angular velocity of the wheel includes the angular velocity of the RS coordinate system. RS

& ω wheel / E = RS ω RS/ E + φ wheel

RS

r$2

(4.3.2)

where φwheel is the angular degree of freedom (scalar) associated with the rear wheel. The value is time dependent and is dictated by the interaction of the vehicle and tire models. The terms in the equations of motion are given by E

ETRT ,qk =

d  ∂TRT  ∂TRT  − dt  ∂q&k  ∂qk

(4.3.3)

Considering just one of the wheels, substituting for the kinetic energy and differentiating gives T

ETRT ,qk

 E d  ∂RS ω wheel / E  ∂RS ω wheel / E  RS  [ J wheel ]RS ω wheel / E =  − ∂qk    dt  ∂q&k T

 ∂RS ω wheel / E  RS  E d RS  +  ω wheel / E   [ J wheel ]  ∂q&k   dt 

80

(4.3.4)

The only degrees of freedom which generate non-zero results for the expression above are the βRS,i and φ wheel. The derivatives of the wheel angular velocities are calculated as follows ∂ RS ω wheel / E ∂ RS ω RS/ E = ∂βRS,i ∂βRS,i

∂RS ω wheel / E ∂RS ω RS/ E = ∂β&RS,i ∂β&RS,i

∂ RS ω wheel / E =0 ∂φwheel E

d RS ω wheel / E = dt

RS

d dt RS d = dt

( (

RS

RS

(4.3.5)

∂ RS ω wheel / E RS = r$2 & ∂φ wheel

& ω RS/ E + φ wheel

RS

)

r$2 +

RS

ω RS/ E ×

(

RS

(4.3.6)

& ω RS/ E + φ wheel

RS

r$2

) (4.3.7)

)

&& ω RS/ E + φ wheel

RS

r$2 +

RS

& ω RS/ E ×φ wheel

RS

r$2 (4.3.8)

E

RS d ∂ RS ω wheel / E E d RS d RS $ = r = ( r$2 )+ 2 & dt ∂φ dt dt wheel

RS

ω RS/ E ×RS r$2

(4.3.9)

= ω RS/ E × r$2 RS

RS

For qk = βRS,i and substituting for the derivatives of

RS

ω wheel / E , using the expressions

above, the result becomes

ETRT ,qk

 E d ∂RS ω T  ∂RS ω T  RS RS/ E RS/ E    [ J wheel ] =  − ∂β &  dt   ∂ β   RS,i   RS, i

(

∂RS ω TRS/ E  RS  RS d   +  [ J wheel ]  dt  ∂β&RS,i 

{

For qk ∈ φleft , φright

(

RS

ω RS/ E

)

RS

RS & r$2 ω RS/ E + φ wheel

 RS & RS r$ + RS ω & r$2  + &φ 2 wheel RS/ E × φwheel 

}the expression simplifies as follows

81

) (4.3.10)

ETRT ,φ =

(

+

RS

RS

(

)

T

ω RS/ E ×RS r$2 [ RS J wheel ]

 RS d r$2T [ RS J wheel ]  dt

(

RS

RS

& RS r$ ω RS/ E + φ 2 right

)

&& RS r$ + ω RS/ E + φ 2 wheel

RS

)

 RS & r$2  ω RS/ E ×φ wheel 

(4.3.11)

4.4 Rear Springs and Dampers The forces exerted by the springs and dampers which connect the rear suspension to the chassis of the vehicle depend on the length and on the rate of change of the length of the springs. The positions of the upper attachment points to the chassis are given by four vectors which are constant with respect to the sprung mass coordinate system SM. Likewise, the lower attachment points are given by four vectors which are constant with respect to the rear unsprung mass coordinate system RS. To determine the terms in the equations of motion which are generated in connection with the spring and damping elements it is necessary to determine the magnitude and direction of the force and then to find the generalized forces. The length of the spring or damper under consideration is E

L i = E r SM / E +

[E CSM ]SM rupper /SM − E r RS/ E − [E C RS ] RS rlower / RS

(4.4.1)

The rate of change of the length is required to determine the force exerted on the system by the dampers. E

& = E r& L i SM / E +

[ C& ]

SM

E

SM

rupper /SM − E r&RS/ E −

[ C& ] E

RS

RS

rlower / RS

(4.4.2)

The forces generated by the springs and by the dampers are determined in the same manner as was done for the front suspension in the previous chapter. For the springs,

82

∆ L = fS − E

FS

(

E

(

E

LS − pS

)(

)

(4.4.3)

)

$ LS = k 0 ∆ L + k1∆ L3 E L S

where E

[L

LS =

E

] 1

TE S

LS

2

(4.4.4)

and for the dampers, E

FD

( L& )= − (c E

D

E R0

)

3 $ L&D + cR1 E L&D E L D

(4.4.5)

where E

d L&D = dt

[L E

T E D

]= 1

LD

2

E

$T EL & L D D

(4.4.6)

The generalized forces are found by determining the virtual work done by the springs and dampers. The total virtual work is

∑

δW =

E

i ∈springs

Fspring, i ⋅δE L i +

∑

E

j∈ dampers

Fdampers, j ⋅δE L j

(4.4.7)

The virtual displacement δE L i is the displacement at the spring or damper caused by an infinitesimal contemporaneous perturbation of the generalized coordinates. The displacement is equivalent to the total derivative of the expression for the length of the spring or damper.

(

δE L i = δE r SM / E +

[C ]

SM

E

SM

rupper /SM − E r RS/ E −

[C ] E

RS

RS

rlower / RS

)

δE L i = (1 0 0)( δ xSM − δx RS ) + ( 0 1 0)( δySM − δ y RS ) + ( 0 0 1)( δzSM − δ z RS ) +

 ∂

SM

∑ ∂β [ C ] E

j

SM

rupper /SM δ βSM, j −

 ∂

E

j

SM, j

83

RS

∑ ∂β [ C ] RS, j

RS

rlower / RSδ βRS, j

(4.4.8)

The generalized forces are equal to the coefficients of the δβSM,j and δβRS, j terms in the expression for the virtual work δW . The generalized forces for a single spring or damper are given below. QRS,xSM = E F ⋅(1 0 0)

QRS,ySM = E F ⋅( 0 1 0)

(4.4.9)

 ∂[E CSM ]SM  rupper /SM QRS,βSM,j = E F ⋅  ∂β   SM, j 

(4.4.10)

QRS,zSM = E F ⋅( 0 0 1)

QRS,x RS = − E F ⋅(1 0 0) QRS,y RS = − E F ⋅(0 1 0)

(4.4.11)

QRS,z RS = − E F ⋅(0 0 1)

∂[E CRS ]RS  QRS,βRS,j = − E F ⋅  ∂β  rlower / RS  RS, j 

(4.4.12)

4.5 Panhard Rod and Trailing Link Constraints The rear suspension is located laterally by a panhard rod and longitudinally by multiple trailing links. These mechanical constraints are represented by an equation which forces the attachment points of the links to remain at a fixed distance from each other. The constraint equation for the panhard rod is derived below. The constraints for the trailing links are identical in form. The vector which lies along the panhard rod is given by the expression E

p = E r SM/ E + [E CSM ]SM rpanhard /SM − E r RS/ E − [E C RS ]RS rpanhard / RS

84

(4.5.1)

The locations of the points of attachment of the panhard rod to the chassis and to the rear suspension are given by the constant vectors

SM

rpanhard /SM and

RS

rpanhard /RS . For a panhard

rod of fixed length L, the following constraint equation can be written.

[ p ⋅ p] − E

E

1 2

L=0

(4.5.2)

The generalized forces associated with the panhard rod constraint are of the form QC ,βi = λC a βi

(4.5.3)

where

aβ i =

1 ∂  E E 12 ∂E p  E  = E p ⋅E p − 2   [ p ⋅ p] − L ⋅ p   [ ]  ∂βi  ∂βi  

(4.5.4)

∂E p = p$ ⋅ ∂βi E

where E p$ is the unit vector along the panhard rod. Evaluating the derivatives of

E

p and

substituting gives QC , xSM = λP E p$ T ⋅(1 0 0)

QC , ySM = λP E p$ T ⋅( 0 1 0)

(4.5.5)

∂[E CSM ] SM  rpanhard /SM QC ,βSM,i = λP E p$ T ⋅  ∂β   SM, i 

(4.5.6)

QC , zSM = λP E p$ T ⋅( 0 0 1)

QC , x RS = − λP E p$ T ⋅(1 0 0) QC , y RS = − λP E p$ T ⋅(0 1 0) QC , z RS = − λP E p$ T ⋅(0 0 1)

85

(4.5.7)

 ∂[ C ] QC ,βRS,i = − λP E p$ T ⋅ E RS   ∂βRS,i 

RS

rpanhard / RS

(4.5.8)

The Lagrange multiplier λP is determined by appending the second time derivative of the constraint equation to the equations of motion. The first and second time derivatives of the constraint equation are calculated below.

{[

] ⇒ [ p ⋅ p] ( d dt

E

}

1

p ⋅E p 2 − L = 0

E

E

−

1 2

E

)

p ⋅E p& = 0

(4.5.9)

⇒ E p ⋅E p& = 0 d E E { p ⋅ p& = 0} dt && + E p&⋅E p& = 0 ⇒ E p ⋅E p

(4.5.10)

The derivatives of the panhard rod vector are computed as follows. E

p = E r SM/ E + [E CSM ]SM rpanhard /SM − E r RS/ E − [E C RS ]RS rpanhard / RS

E

p& = E r&SM/ E +

E

SM

&& = E &r&SM/ E p

E

SM

E

[ C& ] && ] + [C

SM

SM

[ C& ] && ] − [C

rpanhard /SM − E r&RS/ E −

E

RS

rpanhard /SM − E r&&RS/ E

E

RS

RS

rpanhard / RS

RS

rpanhard / RS

4.6 Summary of Results Table 4.1 Kinetic and Potential Energy Terms for the Motion of the Rear Suspension

Generalized Coordinate xSM ySM zSM

Term

mRS x&&RS mRS &y&RS

mRS ( &z&RS + g )

86

(4.5.11)

βRS,i

T

 ∂RS ω RS/ E  RS  d RS    ω RS/ E   ∂β&  [ J RS ]   dt   RS, i +

(

RS

 d ∂RS ω   ∂RS ω RS/ E  T RS/ E    ω RS/ E [RS J RS ] −   dt   &   ∂βRS,i   ∂βRS,i 

)

Table 4.2 Kinetic Energy Terms for the Rotation of the Rear Wheels and Tires

Generalized Coordinate βRS,i

Term  E d ∂RS ω T  ∂RS ω T  RS/ E RS/ E  RS    −   dt  [ J wheel ] & ∂ β  ∂βRS,i  RS, i  

(

RS

RS & r$2 ω RS/ E + φ wheel

 ∂RS ω TRS/ E  RS  E d RS  +  ( ω RS/ E ) + &φ&wheel RS r$2 +  ∂β& [ J wheel ] dt   RS, i 

(

φwheel

+

RS

)

T

ω RS/ E ×RS r$2 [RS J wheel ]

RS

 RS d r$2T [ RS J wheel ]  dt

(

RS

(

RS

& RS r$ ω RS/ E + φ right 2

)

&& RS r$ + ω RS/ E + φ wheel 2

RS

 RS & r$2  ω RS/ E ×φ wheel 

)

RS

 RS & ω RS/ E ×φ r$2  wheel 

Table 4.3 Generalized Forces Associated with a Rear Spring or Damper

Generalized Coordinate xSM, ySM, zSM βSM,i

xRS, yRS, zRS βRS,i

Generalized Force E

E

F ⋅(1 0 0),

E

F ⋅( 0 1 0),

 ∂[E CSM ]SM  rupper /SM F ⋅  ∂β   SM, j 

E

F ⋅( 0 0 1)

− E F ⋅(1 0 0), − E F ⋅(0 1 0), − E F ⋅(0 0 1)

 ∂[E CRS ]RS  rlower / RS − E F ⋅  ∂β   RS, j 

87

)

Table 4.4 Constraint Forces Associated with the Panhard Rod

Generalized Coordinate xSM, ySM, zSM

Generalized Constraint Force

βSM,i

∂[E CSM ] SM  rpanhard /SM λP E p$ T ⋅  ∂β   SM, i 

xRS, yRS, zRS βRS,i

λP E p$ T ⋅(1 0 0), λP E p$ T ⋅( 0 1 0), λP E p$ T ⋅( 0 0 1)

− λP E p$ T ⋅(1 0 0), − λP E p$ T ⋅(0 1 0), − λP E p$ T ⋅(0 0 1)

 ∂  − 12  E  − λ P [E p ⋅E p]  [ E C RS ]  p ⋅  ∂βRS,i  

88

(

RS

rRS/ piv +

RS

 rpanhard / RS   

)

5 Equations of Motion - Steering System 5.1 Introduction The steering system model is tightly integrated with the chassis model and with the front suspension model. Since the input to the steering model is an angular displacement of the steering wheel it is not necessary to consider the inertia of the steering wheel. For this reason, and also to reduce the total number of degrees of freedom, the steering system is modeled as a quasi-static system. If the driver model were to interact with the steering wheel via a prescribed torque it would be better to include the inertia of the wheel in the model. Two types of steering systems are commonly used. The simpler of the two is the rack and pinion type steering system. The second type is based on a four bar linkage.

5.2 Rack and Pinion The rack and pinion steering system is modeled almost trivially. A schematic of the system is shown Figure 5.1. The location of the center of the steering rack relative to the sprung mass coordinate system is given by the vector

SM

R R /SM . The steering rack coordinate

system in this case is simply a displaced version of the sprung mass coordinate system; the

r$1 and r$2 unit vectors are parallel to the corresponding unit vectors of the sprung mass coordinate system. The locations of the tie rod connection points on the rack are given by the relations

89

w

r$1 r$2

R SM

θsw

R R / SM

s$1

SM

s$ 2

Figure 5.1 Schematic of the Rack and Pinion Steering System

SM

w  R LTR/SM = SM R R /SM + G θsw + s$ 2  2

SM

R RTR /SM =

SM

R R /SM

w  + G θsw −  s$ 2  2

(5.2.1)

where G is the gain of the rack and pinion with units of length/angle. The remainder of the steering linkage is represented by a constraint equation which forces the steering knuckle joint on the left and right spindles to remain a fixed distance from the two points given above.

5.3 Four Bar Linkage The four bar linkage type of steering mechanism is considerably more difficult to model than the rack and pinion type mechanism. To facilitate modeling the components of the four bar steering mechanism are assumed to be composed of rigid bodies. The gearbox is modeled as a simple gear reduction without losses. The gearbox output drives the pitman

90

arm which is rigidly attached to the drag link. The pitman arm, idler arm and drag link are modeled as a planar four bar linkage. The position of the four bar linkage is completely specified by the steering wheel angle input. The drag link is connected to the steering knuckles via tie rods which are represented by constraint equations. The plane of the four bar linkage is typically not aligned with the

s$1 − s$ 2 plane of

the sprung mass coordinate system. Given the axis of rotation of the pitman arm (and assuming a parallel axis of rotation for the idler arm) a coordinate system P (for pitman arm) can be defined which contains the four bar linkage in its

p$ 1 − p$ 2 base plane. A

schematic showing the relationship between the sprung mass and pitman arm coordinate systems is shown in Figure 5.2. Note that the view shown in the figure is not necessarily in the s$1 − s$ 2 plane. The orientation of the coordinate system is chosen such that p$ 1 lies in

p$ 1

P

p$ 2

r$PA / IA

P

P P

Pitman Arm

r$PJ / PA

r$IJ /IA Idler Arm

Drag Link

Right Tie Rod Left Tie Rod SM

s$1

r$PA /SM s$3 s$ 2

Figure 5.2: Relationship Between the P and the S Coordinate Systems

91

S

p$ 1 d$ 1 Pitman Arm

p$ 2 d$ 2

D

P

r$LTR / D

D

D

r$RTR/ D

Idler Arm

Figure 5.3: Relationship Between the D and the P Coordinate Systems

the s$1 − s$3 plane and p$ 2 is equal to s$ 2 . Note that this choice only allows the plane of the four bar linkage to be rotated about the s$ 2 axis relative to the sprung mass coordinate system (i.e. only fore-aft tilt of the pitman/idler arm rotational axis is allowed). The vectors P r$PJ / PA and

P

r$IJ / IA represent the pitman arm and the idler arm respectively; the

points PJ and IJ refer to the pitman arm joint and the idler arm joint while points PA and IA refer to the rotational axes of the pitman arm and the idler arm. A second coordinate system D, attached to the drag link, is defined in order to locate the tie rod joints which lie on the drag link. The d$ 2 axis is chosen to lie along the vector r$PJ /IJ (from the idler arm joint on the drag link to the pitman arm joint on the drag link). The unit vector d$ 1 is chosen to be perpendicular to d$ 2 while at the same time remaining in the p$ 1 − p$ 2 plane; it’s direction is chosen to point in roughly the same direction as p$ 1 (in the direction of forward motion). A schematic showing the relationship

92

between the P and D coordinate systems is shown in Figure 5.3. The locations of the tie rod joints on the drag link are given by the vectors

D

r$LTR / D and D r$RTR / D .

The first step in modeling the system is to determine the unit vectors of the P coordinate system which defines the plane of the four bar linkage. The origin of the P coordinate system is on the axis of rotation of the pitman arm and located in the plane of rotation of the pitman arm to drag link joint. Given that the axes of rotation for the pitman and idler arms are assumed to be parallel to the s$1 − s$3 plane the following constraints must be satisfied by the p$ 1 unit vector. SM

p$ 1 ⋅SM p$ 3 = SM p$ 1 ⋅SM a$ axis = 0

(5.3.1)

p$ 1 ⋅SM s$ 2 = 0

(5.3.2)

SM

p$ 1 = 1

The unit vector

SM

(5.3.3)

a$ axis represents the axes of rotation of the pitman arm and idler arm.

Applying the constraints gives the following results.

SM

η 0 − p$ 1 =  

(

SM

a$ axis, x

SM

a$ axis, z

)

1

(

1 SM

a$ axis, x

SM

a$ axis, z

2  2  

)

(5.3.4)

p$ 2 = s$ 2

(5.3.5)

p$ 3 = SM a$ axis

(5.3.6)

SM

SM

  ; η ≡  η  1 + 

To locate the tie rod joints on the drag link it is first necessary to determine the orientation of the drag link and of the idler arm.

93

p$ 1

p$ 1

θP + θP0

SM

SM

r$PA / IA

θI + θI0

r$PJ / PA SM

r$IJ / IA

L DL

Figure 5.4: Schematic of the Four Bar Linkage Steering System

The geometry of the problem is shown in Figure 5.4. The vectors SM

SM

rPJ / PA and

rIJ /IA represent the position of the pitman arm and of the idler arm respectively. The

angles θP0 is the between the p$ 1 unit vector and the pitman arm for zero steer angle (straight ahead steering). θI0 is defined in the same manner but for the idler arm.

θP is

related to the steering wheel angle through a gear reduction. The vectors which represent the pitman arm and the idler arm can be written in terms of the previously defined angles. SM

rPJ / PA = L P [cos(θP + θP0 )SM p$ 1 + sin( θP + θP0 )SM p$ 2 ]

(5.3.7)

rIJ / IA = L I [cos(θI + θI0 )SM p$ 1 + sin( θI + θI0 )SM p$ 2 ]

(5.3.8)

SM

The unknown angle θI , and hence the location of the drag link, can be obtained writing a constraint equation which enforces connectivity between the element of the four bar linkage. SM

rPJ / PA +

SM

rPA /IA −

SM

94

rIJ /IA = L DL

(5.3.9)

Expanding the equation gives the following result:

[ + L [L sin( θ = (L + L + + L [ r

L I L P cos(θP + θP0 )+ I

P

2 P

1 2

2 I

SM

SM

P

+ θP0 )+

P

PA / IA,1

SM

2 rPA / IA,1 +

] ]sin( θ + θ − L )

rPA / IA,1 cos(θI + θI0 )

SM

rPA / IA,2

SM

I

2 rPA / IA,2

cos(θP + θP0 ) +

I0

)

2 DL

SM

(5.3.10)

]

rPA / IA,2 sin( θP + θP0 )

Making the following definitions

[ β= L [ L sin( θ γ= ( L + L + + L [ r

α = L I L P cos(θP + θP0 )+ I

P

P

2 P

1 2

2 I

+ θP0 )+

SM

SM

P

PA /IA,1

SM

SM

2 rPA / IA,1 +

] ]

rPA /IA,1

rPA /IA,2

SM

2 2 rPA /IA,2 − L DL

cos(θP + θP0 ) +

SM

)

(5.3.11)

]

rPA /IA,2 sin( θP + θP0 )

and substituting gives

α cos(θI + θI0 ) + β sin( θI + θI0 ) = γ

(5.3.12)

which can be solved for θI by application of Newton’s method. Using an initial guess of

θI = 0 appears to work well. Once the locations of the pitman arm and of the idler arm are known the vector SM

rPJ / IJ can be determined SM

rPJ / IJ = SM rPJ / PA +

SM

rPA / IA −

SM

rIJ / IA

(5.3.13)

The unit vectors of the D coordinate system can be determined at this point.

SM

d$ 2 =

SM

rPJ /IJ

SM

rPJ /IJ

95

(5.3.14)

SM

(

d$ 1 = cos( φ)SM p$ 1,1 sin( φ) cos( φ)SM p$ 1,3 − φ= arctan    SM

SM

)

p$ 1,1 SM d$ 2,1 + SM p$ 1,3 SM d$ 2,3   SM $  d 

(5.3.15)

2, 2

d$ 3 = SM p$ 3 = SM a axis

(5.3.16)

Finally, the locations of the tie rod joints can be expressed as SM

rLTR/SM =

SM

rRTR / P =

(

(

D

D

rLTR / PJ,1 SM d$ 1

D

rRTR/ PJ,1 SM d$ 1

D

rLTR/ PJ,2 SM d$ 2

rRTR / PJ,2 SM d$ 2

D

D

)

SM

)

SM

rLTR / PJ,3 SM d$ 3 +

rRTR / PJ,3 SM d$ 3 +

rPJ / PA +

SM

rPJ / PA +

SM

rP/SM

(5.3.17)

rP/SM

(5.3.18)

5.4 Tie Rod Constraints Given the locations of the tie rod joints (for either the rack and pinion steering system or the four bar link steering system) found in the previous sections a pair of constraint equations can be written to connect the steering mechanism with the steering knuckles on the spindles. The constraint equation simply expresses the fact that the distance between the ball joint on the steering knuckle and the tie rod joint on the steering mechanism is a constant. That is, E

r KJ /TJ = E r SP/ E +

[C ]

r KJ /SP − E r SM / E −

T E KJ / TJ

r KJ / TJ ] − LTR = 0

SP

[r E

[C ]

SP

E

SM

E

1 2

SM

r TJ /SM

(5.4.1)

(5.4.2)

where TJ indicates the tie rod joint, KJ indicates the steering knuckle joint and LTR is the length of the tie rod. The generalized constraint forces are of the form QSP ,qi = λSP a qi

96

(5.4.3)

where

a qi =

1  ∂  ∂E T E − 12  E T T E [ rKJ / TJ r KJ / TJ ]2 − LTR   = [E rKJ  rKJ / TJ  E r KJ / TJ  r ] /TJ KJ /TJ   ∂qi ∂qi  

{

For qi ∈ xsp , ysp , zsp , xsm , ysm , z sm

∂ E r = ∂xsp KJ / TJ

(1

}the

∂ E ∂qi KJ / TJ

r

∂ E r = ∂y sp KJ / TJ

0 0)

(0

term evaluates to

1 0)

∂ E r = ∂z sp KJ / TJ

(0

0 1)

∂ E ∂ E ∂ E r KJ / TJ = − (1 0 0) r KJ / TJ = − ( 0 1 0) r = − ( 0 0 1) ∂xsm ∂ysm ∂z sm KJ / TJ

{ }

∂ E ∂qi KJ / TJ

For qi ∈ βSP,i the

 SP  ∂ ∂ E  r KJ /SP r KJ / TJ =  [ E C SP ] ∂βSP,i   ∂βSP,i

{

}

(5.4.5)

term evaluates to

r

For qi ∈ βSM, i the

(5.4.4)

∂ E ∂qi KJ / TJ

(5.4.6)

term evaluates to

r

 SM  ∂ ∂ E  r TJ /SM r KJ /TJ = −  [ E C SM ] ∂βSM, i   ∂βSM, i

(5.4.7)

In addition to the generalized forces the second time derivative of the constraint equation is required in order to solve for the Lagrange multiplier. The first derivative is calculated as follows

d dt ⇒

{[ r E

T E KJ / TJ

[r E

T E KJ / TJ

}

r KJ / TJ ] − LTR = 0 1 2

r KJ / TJ ] −

1 2

E

T E & rKJ / TJ r KJ / TJ = 0

T E & ⇒ E rKJ / TJ r KJ / TJ = 0

and the second derivative is

97

(5.4.8)

d E T E { rKJ / TJ r&KJ / TJ = 0} dt T E E T E & && ⇒ E r&KJ / TJ r KJ / TJ + rKJ / TJ r KJ / TJ = 0

(5.4.9)

The first and second time derivatives of the position vector are simply E

E

r&KJ /TJ = E r&SP/ E +

[ C& ]

r KJ /SP − E r&SM / E −

[ C& ]

r TJ /SM

(5.4.10)

r&&KJ /TJ = E r&&SP/ E +

[ C&& ]

r KJ /SP − E &r&SM / E −

[ C&& ]

r TJ /SM

(5.4.11)

SP

E

SP

SP

E

SP

98

SM

E

SM

SM

E

SM

6 Road Model 6.1 Introduction The road model is a critical component of the vehicle simulation. If the complete simulation is to accurately represent reality the road model must accurately represent the terrain. In addition to modeling the large scale features of a particular road course the road model must also be capable of representing the small scale features such as bumps and other road surface irregularities. On the other hand, it is undesirable for the road model to be so detailed that massive amounts of data are required to generate a complete road course. It is necessary to balance the competing demands of accuracy and a minimal data set to obtain a useful road model. The model described below is an effective compromise in attaining these opposing goals. The model utilizes parametric polynomial equations in three dimensions to describe the path followed by the centerline or the road. Polynomials of various orders are used depending on the nature of the boundary conditions on a particular segment of the road. In addition to specifying the endpoints (and the tangent vectors at the endpoints for the higher order polynomials) it is also necessary to provide the angle of the road surface normal with respect to the vertical. This surface normal angle specifies the lateral tilt of the road. It is assumed that the road is flat in the lateral direction (no crown or concavity of the road surface). 99

The modeling methodology laid out in the preceding paragraph is sufficient for describing the large scale features of the roadway and perhaps even some of the larger bumps on the road surface. The description of small scale features requires an addition to the model. The texture associated with the road surface can be simulated by a filtered white noise function. A relatively small number of parameters are required to describe the average amplitude of the road surface irregularities as a function of their size. Prior to actually running the simulation, the irregularities can be generated for each road segment and tabulated for quick lookup. The number of entries in the lookup table is determined by the size of the smallest irregularity being modeled.

6.2 Road Surface Coordinate System As discussed above the road is modeled by specifying a parametric function which gives the location of the road centerline and by specifying an angle which gives the tilt of the road surface from the vertical. It is useful to determine a set unit vectors (in terms of the inertial coordinate system) which defines a road surface coordinate system. This is most easily done by utilizing an intermediate coordinate system T. The intermediate coordinate system is defined such that the x-axis unit vector the path. The y-axis unit vector

E

E

t$ x lies along the tangent to

t$ y is perpendicular to the x-axis vector and is parallel to

the x-y plane of the inertial system. E

t$ x = E rP/′E ( s ) / E rP/′E ( s )

100

(6.2.1)

where and

E

E

rP/E ( s) is the parametric equation specifying the location of the road center line

rP/E ′( s) is the tangent to the path (obtained by differentiating

E

rP/E ( s) with respect

to s). The y-axis unit vector is of the form E

(

$t y = t y1

)

t y2

0

(6.2.2)

The orthogonality condition with tx dictates that t$

E TE x

$t y = t x1t y1 + t x 2 t y 2 = 0

(6.2.3)

and the normalization condition requires that t y21 + t y22 = 1

(6.2.4)

Solving for the components of ty gives

t y1 = ±

t x22 t x21 + t x22

t  t y 2 = −  x1 t y1 t x 2 

(6.2.5)

In the event that tx2 is much less than tx1 (or even if it is zero) then the alternative solution below may be used.

t y1

t  = −  x 2 t y 2  t x1 

t y2

t x21 =± 2 t x1 + t x22

(6.2.6)

The sign of the solution is chosen such that the z-axis of the intermediate coordinate system points up (i.e. in roughly the same direction as the z-axis of the inertial coordinate system). E

t$ z = E t$ x ×E t$ y

tz3 > 0

(6.2.7)

The road surface coordinate system is obtained from the intermediate coordinate system via a simple rotation which accounts for the lateral tilt of the road surface:

101

E

n$ x = E t$ x

E

n$ y = cos(θ)E $t y − sin( θ)E t$ z

E

n$ z = sin( θ)E t$ y + cos(θ)E $t z

(6.2.8)

Note that θ may be a function of the free parameter.

6.3 Location of Tire to Road Contact Point In order to model the interaction of the road surface and the vehicle model it is necessary to determine the location of the point of contact between the tire and the road. The contact point is known to lie in the plane of the road and is assumed to lie in the plane of the wheel. The intersection of these planes is a line. The point on this line which minimizes the distance to the center of the wheel is taken to be the contact point. These conditions form a constrained minimization problem which can be solved using the method of Lagrange multipliers.

f ≡ min E rCP/E − E rWC/E

(6.3.1)

c1 ≡ ( E rCP/ E − E rWC/ E ) E n$ W = 0

(6.3.2)

c2 ≡E rCP/ E = α E n$ y ( s )+ E rP/ E ( s )

(6.3.3)

T

The variables are defined as follows:

E

rCP/E is the unknown location of the contact point,

E

rWC/E is the location of the wheel center with respect to the inertial coordinate system,

E

n$ W is the normal to the plane of the wheel, α is an unknown constant which is equal to

the distance of the contact point from the road center line, s is the unknown parameter giving the location along the path, and

E

n$ y ( s ) is the normalized y-axis vector of the road

surface coordinate system. Rather than minimizing the expression

102

E

rCP/E − E rWC/ E , which

(

contains a square root, it is simpler to minimize

rCP/E − E rWC/E )

T

E

(

E

rCP /E − E rWC /E ). The

system of equations can be further simplified by substituting the last equation into the first two equations; for the moment though, the first equation is left in terms

E

rCP/ E to allow

easier manipulation.

f ≡ min ( E rCP/ E − E rWC / E )

(

T E

rCP/ E − E rWC/ E )

(

c1 ≡ α E n$ y ( s )+ E rP/ E ( s )− E rWC / E

)

TE

n$ W = 0

(6.3.4)

Applying the method of Lagrange multipliers gives a system of three of equations

∇ f + λ∇ c1 = 0

(6.3.5)

c1 = 0

(6.3.6)

This system of equations has two unknowns: α and s. Thus, the del operator is defined as

∇ = ( ∂∂α

)

∂ T ∂s

(6.3.7)

Evaluating the derivatives in the preceding equations gives

∂ ∂α

c1 =

∂ ∂α

∂ ∂s

∂ ∂α

f = 2( E rCP/ E − E rWC/ E )

(

∂E ∂α CP/ E

r

)

(6.3.8)

∂ ∂s

f = 2( E rCP/ E − E rWC/ E )

(

∂E ∂s CP/ E

)

(6.3.9)

(α

T

T

E

n$ y ( s )+ E rP/E ( s )− E rWC/E

(

r

)

TE

n$ W = E n$ TW E n$ y ( s )

c1 ≡ ∂∂s α E n$ y ( s )+ E rP/ E ( s )− E rWC/ E

[(

= α

∂E ∂s

) (

n$ y ( s ) +

∂E P/ E ∂s

r

)]

T

(s)

)

TE

E

(6.3.10)

n$ W (6.3.11)

n$ W

where ∂E CP/ E ∂α

r

=

∂ ∂α

(α

E

)

n$ y ( s )+ E rP/E ( s ) = E n$ y ( s )

103

(6.3.12)

∂E CP/ E ∂s

r

The terms

∂E ∂s

n$ y ( s ) and

=α

∂E ∂s P/ E

r

(

∂E ∂s

) (

n$ y ( s ) +

)

∂E P/E ∂s

r ( s)

(6.3.13)

( s ) depend on the form of the parametric function being

used to model the road segment. Expressions for these terms are derived in sections below. Assembling terms gives the final system of equations which can be solved to determine the Lagrange multiplier, α, and s. The location of the contact point is determined trivially from α and s.

[(

)

f 0 ≡ 2 α E n$ y ( s )+ E rP/ E ( s )− E rWC/ E + λE n$ W

[(

)

f1 ≡ 2 α E n$ y ( s )+ E rP/ E ( s )− E rWC / E + λE n$ W

] (α( T

[

]

T

E

n$ y ( s ) = 0

) (

∂E ∂s

n$ y ( s ) +

]

n$ W = 0

f 2 ≡ α E n$ y ( s )+ E rP/ E ( s )− E rWC/ E

TE

(6.3.14)

∂E P/ E ∂s

r

))

(s ) = 0

(6.3.15)

(6.3.16)

The three equations above are typically nonlinear for anything but the simplest of road segment geometries. To obtain the solution it is necessary to apply a multidimensional form of Newton’s method or a similar algorithm. Most of the efficient algorithms require evaluation of the Jacobian. The terms of the Jacobian are computed below:

∂f 0 = 2E n$ Ty E n$ y ∂α

(6.3.17)

 ∂En$ y ∂ErP/ E  ∂ f1  + 2 α E n$ y + E rP/ E − E rWC/ E + λE n$ W = 2 E n$ Ty α + ∂α ∂s   ∂s

[(

)

∂f 2 E T E = n$ W n$ y ∂α

104

]

T

∂En$ y ∂s

(6.3.18)

(6.3.19)

∂f 0 E T E = n$ W n$ y ∂λ

(6.3.20)

E ∂f1 E T  ∂ n$ y ∂ErP/ E   = n$ W α + ∂λ ∂s   ∂s

(6.3.21)

∂f 2 =0 ∂λ

(6.3.22)

 ∂En$ y ∂ErP/ E  ∂f 0  = 2α + ∂s ∂s   ∂s

T E

[(

)

n$ y + 2 α E n$ y + E rP/ E − E rWC/ E + λE n$ W

]

T

∂En$ y ∂s

(6.3.23)

T

 ∂En$ y ∂ErP/ E   ∂En$ y ∂ErP/ E  ∂ f1  α  = 2α + + ∂s ∂s   ∂s ∂s   ∂s

[(

)

+ 2 α E n$ y + E rP/ E − E rWC/ E + λE n$ W

]

T

(6.3.24)

 ∂2 E n$ y ∂2 E rP/ E  α  + 2 ∂s 2   ∂s

E ∂f 2 E T  ∂ n$ y ∂ErP/ E   = n$ W α + ∂s ∂s   ∂s

The direction vector

E

(6.3.25)

n$ y ( s ) has not yet been determined. The direction vector for the x-

axis of the road coordinate system is given by (recalling the results from the preceding section) E

n$ x ( s )= E t$ x ( s )= E rP/′E ( s ) / E rP/′E ( s )

where the prime indicates differentiation with respect to s. Let perpendicular to

E

(6.3.26) E

t$ y ( s ) be a vector

$t x ( s ) and oriented such that it lies in the x-y plane of the inertial

coordinate system (again, following the development in the first section). If form

105

E

$t x ( s ) has the

E

T t$ x ( s ) = (t x1 ( s ) t x 2 ( s ) t x 3 ( s ))

(6.3.27)

then E

(

t$ y ( s ) = t y1 ( s ) t y 2 ( s ) 0 =

1 t x21 ( s ) + t x22 ( s )

)

(−

T

(6.3.28)

t x 2 ( s ) t x 1 ( s ) 0)

T

and E

T t$ z ( s ) = (t z1 ( s ) t z 2 ( s ) t z 3 ( s ))

(

= E $t x ( s )×E t$ y ( s ) = − t x 3t y 2

t x 3 t y1

t x1t y 2 − t x 2 t y1

)

T

(6.3.29)

where the tx,i depend on the form of the parametric equation used for a particular road segment. The tx,i are derived for several types of road segments in the following section. The particular orientation of E

E

t$ y ( s ) is irrelevant (to the left of

E

$t x ( s ) or to the right of

t$ x ( s ) ); the only result will be a change of sign of α in the solution. Applying the rotation

to account for the tilt of the road gives

E

 t y1 cos(θ ) + t z1 sin( θ )    n$ y ( s ) = t y 2 cos(θ ) + t z 2 sin( θ )   t z 3 sin( θ )  

(6.3.30)

where θ is typically a function of s. The first and second derivatives of

E

n$ y ( s ) are also required:

 ∂∂t ys1 cos(θ ) − t y1 sin( θ ) ∂θ  + ∂∂t zs1 sin( θ ) + t z1 cos(θ ) ∂θ ∂s ∂s   ∂E ∂t y 2 ∂t z 2 ∂θ  n$ y ( s ) =  ∂s cos(θ ) − t y 2 sin( θ ) ∂θ + sin( θ ) + t cos( θ ) z2 ∂s ∂s ∂s ∂s   ∂t z 3 ∂θ   sin( θ ) + t z 3 cos(θ ) ∂s ∂s  

106

(6.3.31)

( ) ( ) ( ) ( ) ( )

 ∂2 t y1 cos(θ ) − 2 ∂t y 1 sin(θ ) ∂θ − t cos(θ ) ∂θ 2 − t sin( θ ) ∂2 θ  y1 y1 ∂s ∂s ∂s  ∂s 2 ∂s 2    ∂2 t ∂t ∂θ ∂θ 2 ∂2 θ  + ∂s 2z 1 sin( θ ) + 2 ∂sz 1 cos(θ ) ∂s − t z1 sin( θ ) ∂s + t z1 cos(θ ) ∂s 2    ∂2 t ∂t y 2 ∂θ ∂θ 2 ∂2 θ y2 sin( θ ) − t 2   2 cos(θ ) − 2 ∂s sin( θ ) ∂s − t y 2 cos(θ ) ∂s 2 y2 ∂s ∂s ∂ E $   n ( ) = s y ∂s 2 + ∂2 t 2z 2 sin(θ ) + 2 ∂t z 2 cos(θ ) ∂θ − t sin( θ ) ∂θ 2 + t cos(θ ) ∂2θ2  z2 z2 ∂s ∂s ∂s ∂s   ∂s 2 2 2 t t ∂ ∂ ∂θ ∂ θ ∂ z 3 z 3  sin( θ ) + 2 ∂s cos(θ ) ∂s − t z 3 sin( θ ) ∂s + t z 3 cos(θ ) ∂s θ2    ∂s 2      

(6.3.32)

where

[ ] [ ] ∂θ = [3θ ′+ 3θ ′ − 6(θ − θ )]s + [6(θ − θ ) − 4θ ′+ 2θ ′ ]s + θ ′ ∂s

θ( s ) = θ0′+ θ1′ − 2(θ1 − θ0 ) s 3 + 3(θ1 − θ0 ) − 2θ0′+ θ1′s 2 + θ0′+ s θ0 (6.3.33)

2

0

1

1

0

1

0

0

1

0

∂2 θ = 6 θ0′+ θ1′ − 2(θ1 − θ0 ) s + 6(θ1 − θ0 )− 4θ0′+ 2θ1′ ∂s 2

[

The derivatives of

E

]

t$ x , E t$ y and

E

t$ z are also required. The derivatives of

E

t$ x depend on

the functional form of the particular type of road segment begin considered. The derivatives of E t$ y and E t$ z have a specific functional dependence on the derivatives of

E

t$ x

and can be determined here.

(

− ∂E $ ty = ∂s

∂t x 2 ∂s

)

∂t x1 ∂s

t + t 2 x1

0 2 x2

T

−

t x1 ∂∂t xs1 + t x 2

[t

2 x1

+ t

107

2 x2

∂t x 2 ∂s

] 3

2

(−

tx2

t x1

0)

T

(6.3.34)

(

− ∂2 E $ t = y ∂s 2

∂2 t x 2

∂2 t x 1

∂s 2

∂s 2

)

T

0

− 2

t x21 + t x22

t x1 ∂∂t xs1 + t x 2

[t

2 x1

+ t x22

∂t x 2 ∂s

] 3

2

+ ( ) ( )+t − (− t + t t [ ] ( t + t ) ( − t t 0) + 3 [t + t ] 2

t x1 ∂∂st2x 1 +

∂2 t x 2 x 2 ∂s 2

∂t x 1 2 ∂s 2 x1

2 x2

∂t x 2 2 ∂s

2

∂t x 2 2 x 2 ∂s

∂t x 1 x 1 ∂s

2 x1

3

2 x2

5

(− x2

∂t x 2 ∂s

∂t x 1 ∂s

t x1

0)

0

)

T

T

(6.3.35)

T

x2

2

x1

∂t  − ∂∂t xs3 t y 2 − t x 3 ∂ys2  ∂E $ ∂t ∂t x 3 t z ( s) =  t + t x 3 ∂sy 1 ∂s y1 ∂s   ∂t x 1 t y 2 + t x1 ∂t y 2 − ∂t x 2 t y1 − t x 2 ∂s ∂s  ∂s

   ∂t y 1   ∂s 

∂t ∂t ∂2 t x 3  t y 2 − 2 ∂∂t xs3 ∂ys2 − t x 3 ∂s y2 2 − ∂s 2  ∂2 E $  ∂2 t y 1 ∂2 t x 3 ∂t x 3 ∂t y 1 t tz = + 2 t y1 + 2 ∂s x 3 ∂ s 2 ∂s ∂s 2  ∂s 2 2 2 ∂ ∂ t t ∂ ∂ ∂ ∂t x 2 t t t y2  x1 t + 2 x1 y 2 + t x2 x 1 ∂s 2 − ∂s 2 t y 1 − 2 ∂s ∂s ∂s  ∂s 2 y 2

(6.3.36)

2

∂t y 1 ∂s

    ∂2 t y1  − t x 2 ∂s 2 

(6.3.37)

6.4 Velocity of the Tire to Road Contact Point The velocity of the tire to road contact point is required by the tire model. The velocity can be found by differentiating the system of equations developed in the preceding section.

d f0 ⇒ dt

[2(

E

)

& n$ + λ n$& r&CP/ E − r&WC / E + λ W W E

d f1 ⇒ dt

[2( r& + [2( r E

CP/ E

E

E

)

]

T

E

[(

&E n$ + λE n&$ − E r&WC / E + λ W W

E

CP/ E

)

− E rWC/ E + λE n$ W

] (α&( T

)

n$ y + 2 rCP/ E − rWC / E + λ n$ W

] (α( T

∂E ∂s

E

∂E ∂s

E

) (

n$ y + α

108

) (

n$ y +

∂2 E ∂s 2

∂E P/ E ∂s

r

E

))

) (

n$ y s&+

]( T

∂E ∂s

)

(6.4.1)

n$ y s& = 0

(6.4.2)

∂2 E P/ E ∂s 2

r

)s&) = 0

d f2 ⇒ dt

[ r& E

CP/ E

− E r&WC/ E

]

TE

n$ W +

[r E

CP/ E

− E rWC / E

]

TE

$W = 0 n&

(6.4.3)

In the preceding section the position of the contact point was specified in terms of a point on the road centerline and the y-axis unit vector road coordinate system. Differentiating that expression gives E

r&CP/ E = α& E n$ y + α

(

∂E ∂s

) (

n$ y s&+

∂E P/ E ∂s

r

)s&

(6.4.4)

The derivatives of the contact point velocity with respect to α& and s& are ∂E CP/ E ∂s&

r&

=α

(

∂E ∂s

∂E CP/ E ∂α&

r&

) (

n$ y +

∂E P/ E ∂s

r

)

(6.4.5)

= E n$ y

(6.4.6)

Once the position of the contact point has been determined using the procedure from the preceding section the only remaining unknowns in the differentiated system of equations

& . Unlike the problem of the preceding section the unknowns appear are α& , s& and λ linearly in the equations and standard linear systems techniques can be applied to determine the solution. A little manipulation of the three equations gives    2 E n$ Ty E n$ y  E T E  n$ W n$ y  T E  2α ∂ E n$ + 2 ∂ E r n$ y P/ E y ∂s ∂s  T  E E E $ + 2 − + r r n λ CP/ E WC/ E W  

( ( [(

))

) (

)

](

T

         E ∂ n$ y  ∂s 

)

109

α&  & E E& λ = 2 r&WC / E − λ n$ W  s&  

[

]

T

E

n$ y

(6.4.7)

((

))

) (

 2 E n$ Ty α ∂∂s E n$ y + ∂∂s E rP/ E   T  + 2 E rCP/ E − E rWC/ E + λE n$ W ∂∂s E n$ y  E T n$ W α ∂∂s E n$ y + ∂∂s E rP/ E   T 2 α ∂∂s E n$ y + ∂∂s E rP/ E α ∂∂s E n$ y + ∂∂s E rP/ E  T 2 2 E + 2 E r − E r $ + λ n α ∂∂s 2 E n$ y + ∂∂s 2 E rP/ E CP/ E WC / E W  

[(

(( [(

)

((

](

)

) ( )) )) ( ( ) ( ))

) (

]( (

)

[

$W = 2 E r&WC/ E − λE n&

] (α( T

∂E ∂s

) (

n$ y +

)(

∂E P/ E ∂s

r

T

           

α&  & λ  s&  

(6.4.8)

))

))

⋅

   α 

[(

n$ Ty E n$ W 0 + ∂∂s E rP/ E

E

∂E ∂s

n$ y

) (

T

   T E n$ W  

)]

α&  & E T E λ= r&WC/ E n$ W −  s&  

[r E

CP/ E

− E rWC/ E

]

TE

n&$ W

(6.4.9)

6.5 Vehicle Position and Heading Angle In order for the vehicle to be driven along the road it is necessary to implement some sort of steering control. Many of the steering control algorithms in use today require feedback of the vehicle location error and the vehicle heading error (or sometimes the vehicle orientation error) with respect to the road. Expressions for the vehicle position and the vehicle heading are developed below. The position of the vehicle relative to the road is a somewhat difficult to quantify in that there are numerous points on the vehicle which may be used as a reference. To simplify matters a single point on the sprung mass can be chosen as a reference point. For most situations it seems reasonable to chose a point in the longitudinal plane of symmetry of the vehicle, somewhere in the vicinity of the center of gravity. Varying the fore/aft

110

location of the point will likely have an effect on the stability and or sensitivity of the control algorithm and so the point should be chosen carefully. The chosen point will most likely be a significant distance above the road surface. Since the distance of the reference point above the road surface is of little concern in a steering control algorithm it is desirable to determine only the lateral components of the position. This can be done by projecting the reference point onto the road surface. Finding the location of the projected point is a good deal easier than finding the contact point between the road and the tire. A set of equations can be written expressing the desired constraints on the projected point: E

rPRJ / E = α E n$ y + E rP/ E

(6.5.10)

(

E

rREF/ E − E rPRJ / E

)

n$ x = 0

(6.5.11)

(

E

rREF/ E − E rPRJ / E

)

n$ y = 0

(6.5.12)

TE

TE

The first equation simply specifies that the projected point lies on the road surface while the second and third equations specify that the vector between the reference point and the projected point is normal to the surface. Substituting the first equation into the second and eliminating the term which contains the product of orthogonal unit vectors gives

(

E

rREF/ E − E rP/ E

)

TE

n$ x = 0

(6.5.13)

This equation no longer contains the unknown α although for most types of road segment models it is strongly nonlinear in s. Once s has been determined the third equation can be used to determine α. α=

(

E

rREF/ E − E rP/ E

111

)

TE

n$ y

(6.5.14)

The velocity of the projected point is determined by differentiating the preceding equations with respect to time: E

(

E

v REF/ E −

α& =

(

E

r&PRJ / E = α& E n$ y + α

(

r

)s&)

(

r

∂E P/ E ∂s

v REF/ E −

T

∂E P/ E ∂s

E

(

(

n$ x +

)s&)

T

E

∂E ∂s

E

) (

n$ y s&+

∂E P/ E ∂s

)s&

)(

∂E ∂s

r

rREF/ E − E rP/ E

n$ y +

(

E

T

rREF/ E − E rP/ E

(6.5.15)

)

n$ x s& = 0

)( T

∂E ∂s

)

n$ y s&

(6.5.16)

(6.5.17)

Rearranging the second equation to isolate s& gives

s& =

E

[(

)

E ∂E T P/ E ∂s

r

n$ x −

v TREF/ E E n$ x

(

E

rREF/ E − E rP/ E

)( T

∂E ∂s

n$ x

)]

(6.5.18)

The second and third equations can be solved trivially once α and s have been determined. The velocity of the projected point is then found using the first equation. It is also useful to determine the error in the heading angle. The heading angle error is determined by projecting the longitudinal axis of the sprung mass coordinate system onto the road surface finding the angle between it and the tangent vector to the road. E

s$ x = [E CSM ](1 0 0)

T

E E

h$ x =

E

( −(

s$ x −

E

s$ x

E

) s$ s$ ) s$

n$ Tz E s$ x n$

E

x

TE E z x x

(

ψ = acos E n$ xT E h$ x

(6.5.19)

)

6.6 Road Segment Models A number of road segment models are formulated in the following sections. The terms which are required in the equations from the preceding sections are calculated for

112

each type of segment. All of the segments derived below are based on polynomial functions of various orders to match different types of boundary conditions. It is possible to base road segments on other types of functions. Note that for all road segment models the free parameter (denoted by s) is assumed to range from zero to one. Also note that the equation expressing the tilt of the road surface is linear for all road segment models and takes the form θ( s ) = ( θ1 − θ0 )s + θ0

(6.6.1)

where θ0 is the tilt of the road at the beginning of the segment ( s=0) and θ1 is the tilt of the road at the end of the road segment (s=1).

Linear Polynomial Road Segment The linear polynomial road segment is based on a parametric equation of the form E

rP/ E ( s ) = a s + b

(6.6.2)

There are two unknown vectors which can be determined via a pair of boundary conditions. The boundary conditions are as follows: E

rP/ E ( 0) = r0

E

rP/ E (1) = r1

(6.6.3)

b = r0

(6.6.4)

Solving for the unknowns gives a = r1 − r0

The tangent vector is obtained by differentiation with respect to s E

r − r t$ x ( s ) = 1 0 r1 − r0

(6.6.5)

The overall path length along the segment is determined trivially in this case L = r1 − r0 113

(6.6.6)

Quadratic Polynomial Road Segment The quadratic polynomial road segment is based on a parametric equation of the form E

rP/ E ( s ) = a s 2 + b s + c

(6.6.7)

The first two boundary conditions are identical to those for the linear road segment. The remaining boundary condition allows matching the tangent vector at one of the endpoints. The choice of which endpoint to use leads to two possible cases. For the first case the boundary conditions are E

rP/ E ( 0) = r0

E

rP/′E (0) = t 0

rP/ E (1) = r1

(6.6.8)

E

rP/′E (1) = t1

(6.6.9)

c = r0

( case 1)

(6.6.10)

c = r0

(6.6.11)

E

and for the second case the boundary conditions are E

rP/ E ( 0) = r0

E

rP/ E (1) = r1

Solving for the unknowns gives a = r1 − r0 − t 0 a = t1 − (r1 − r0 )

b = t0

b = 2(r1 − r0 ) − t 0

( case 2 )

where the tangent vector is obtained by differentiation with respect to s E

$t x ( s ) = 2a s + b 2a s + b

(6.6.12)

The derivatives of E $t x ( s ) can be shown to be

[(

) ]

2 a − E t$ xT a E t$ x ∂E $ t x (s) = ∂s 2a s + b

(

2 ∂2 E $ t x (s) = − 2 ∂s

)

(

t$ a E t$ x + 4 E t$ xT a

∂E T x ∂s

2a s + b

(6.6.13)

)

t$

∂E x ∂s

(6.6.14)

The overall path length along the segment is determined by evaluating the following integral

114

s

L( s ) = ∫

E

T E rP/′ ′E ( s ) ds E ( s ) rP/

(6.6.15)

0

While a closed form solution to this integral may be possible, numerical integration appears to be preferable.

Cubic Polynomial Road Segment The cubic polynomial road segment is based on a parametric equation of the form E

rP/ E ( s ) = a s 3 + b s 2 + c s + d

(6.6.16)

A total of four boundary conditions are required to determine the coefficients. The first two boundary conditions force the function to pass through the starting and ending points and are identical to those for the linear road segment and the quadratic road segment. The remaining two boundary conditions force the tangent vector at the end points to match a specified value. The boundary conditions are E

rP/ E ( 0) = r0

E

rP/′E ( 0) = t 0

E

rP/ E (1) = r1

E

rP/′E (1) = t1

(6.6.17)

Solving for the unknowns gives a = t1 + t 0 − 2(r1 − r0 )

b = − t1 − 2t 0 + 3(r1 − r0 )

c = t0

d = r0

(6.6.18)

where the tangent vector is obtained by differentiation with respect to s

E

2 $t x ( s ) = 3a s + 2b s + c 3a s 2 + 2b s + c

(6.6.19)

The derivatives of E $t x ( s ) can be shown to be

(

)

(6as + 2b)− E $t Tx (6as + 2b) E t$ x ∂E $ t x (s) = ∂s 3a s 2 + 2b s + c

115

(6.6.20)

6a − ∂2 E $ t x (s) = 2 ∂s

(

)

(

)(

t$ (6as + 2b)+ 6 E $t Tx a E t$ x − 2 E $t Tx (6as + 2b)

∂E T x ∂s

t$

∂E x ∂s

3a s + 2b s + c 2

)

(6.6.21)

The overall path length along the segment is determined by evaluating the following integral in the same manner as was done for the quadratic road segment. s

L( s ) = ∫

E

T E rP/′ ′E ( s ) ds E ( s ) rP/

(6.6.22)

0

There is no readily available closed form solution for the resulting integral. Numerical integration appears to be the only viable means to obtain a solution.

116

7 Equations of Motion - Tire Model 7.1 Introduction The tire model plays an important role in the performance of the vehicle model. With the exception of aerodynamic forces which are only relevant at higher speeds, all of the interactions of the vehicle with the external environment occur through the tire. Accurate modeling of the tire is crucial to an accurate vehicle simulation. In keeping with current modeling practices the model for horizontal force generation is decoupled from model for vertical force generation with the exception of the dependence of horizontal force on normal load. The tire model is thus divided into two components. The first component models the vertical support characteristics of the tire utilizing a simple nonlinear spring. Damping in the tire is neglected due to its small magnitude relative to the damping in the shock absorber. The second model component handles the generation of lateral and longitudinal forces using the normal force determined by the vertical model.

7.2 Coordinate Systems There are two coordinate systems which are utilized in tire modeling. They are depicted in Figure 7.1. The wheel coordinate system W is chosen such that the plane of symmetry of the wheel lies in the x-z plane. The x-axis lies in the ground plane. The tire

117

model coordinate system T is defined

$z w

$t z

by choosing $t z normal to the road Wheel Plane

surface at the point of contact and choosing $t x to point in the direction T, W

of the wheel heading ( not necessarily $x w t$ x

the direction of the wheel velocity). The forces and moments which are

$t y

Figure 7.1: The Tire Model Coordinate Systems

determined by the tire model are referred to the T coordinate system. Before proceeding with the actual tire modeling it is necessary to determine the unit vectors of the W and T coordinate systems. The T coordinate system is closely related to the road surface coordinate system described in the chapter on road modeling. The difference is a single rotation about the road surface normal ( $t z axis) to bring the $t x axis into alignment with the wheel heading vector. Given a set of unit vectors r$i which define the road surface coordinate system and the unit vectors for the spindle coordinate system p$ i the T coordinate system unit vectors are determined as follows: E

E

t$ x =

$t z = E r$z

( −(

(7.2.1)

) $t $t ) $t

E

p$ x −

E

p$ Tx E $t z

E

p$ x

E

p$

E

t$ y = E t$ z ×E t$ x

E

z

TE E x z z

118

(7.2.2)

(7.2.3)

It is assumed in the above equations that p$ x points in a direction which is roughly in the direction of travel. The equations above have a singularity when p$ x coincides with $t z . The W coordinate system is trivially derived from the unit vectors of the spindle coordinate system and the unit vectors of the T coordinate system:

E

E

$ x = E t$ x w

(7.2.4)

E

$ y = E p$ y w

(7.2.5)

$ z =Ew $ x ×E w $y w

(7.2.6)

It is also useful to calculate the camber angle of the wheel for use in the horizontal tire force calculations. The camber angle is simply the angle between the

$ z unit vector w

and the $t x − $t z plane. The magnitude can be found by taking the inverse cosine of the

$ z unit vector lies in the $t y − $t z plane. The sign is determined by inner product since the w $ z along $t y . Camber is taken to be positive for positive looking at the component of w rotations about the $t x axis (using the right hand rule).

(

E TE  $ z ≥ 0, - cos-1 E $t zT E w $z  t$ y w γ=  E T E -1 E $ T E $ $ $z tz w   t y w z < 0, cos

(

)

)

(7.2.7)

7.3 The Magic Formula Tire Model The tire itself is modeled using a variation on the popular Magic Formula Tyre Model [Bakker, Pacejka]. The formula is based on a function whose behavior approximates the shape of the curves obtained from experimental measurements on tires. It’s parameters are determined so as to fit the curve to a particular set of experimental

119

data. The function has the added benefit that the coefficients describing the shape of the curve have simple interpretations. The tire model can be divided into two sub-models: one for the vertical (support) force and one for the horizontal (tractive) forces. The only coupling between the two sub-models is the dependence of the tractive forces on the normal load. The support force has no dependence on the tractive force.

Support Forces The support force generated by the tire can be modeled by a nonlinear spring which is placed between the hub of the wheel and the point of contact between the tire and the road surface. Damping effects are similarly modeled using a nonlinear viscous type damping element. The vector from the tire-to-ground contact point to the center of the wheel must be found in order to determine the length of the spring. The time derivative of the same vector must also be determined in order to find the magnitude of the damping force. The calculations for determining the position and velocity of the point of contact between the tire and the road are derived in the chapter on road modeling. The direction of the force is assumed to be along the surface normal and camber effects are ignored (although they appear in the tractive force model). Once the location of the contact point

E

rCP/ E is known the generalized forces

associated with the tire can be determined. The vector between the contact point and the wheel center is given by E

rWC/ CP = E rWC / E − E rCP/ E = Re E n$ WC /CP

120

(7.3.1)

where

E

n$ WC/ CP is the unit vector along

(equal to the magnitude of

E

E

rWC/ CP and Re is the effective radius of the tire

rWC/ CP ). The time derivative of this vector is simply r&WC/ CP = E r&WC/ E − E r&CP/ E

E

(7.3.1)

The rate of change in the length of the vector is determined as follows: d R&e = dt

[r E

T E WC /CP WC /CP

r

]= 1

2

E

n$ TWC /CP E r&WC/ CP

(7.3.1)

Thus, the force exerted by the spring on the wheel center is E

[

Ftire,z = k 0 ( R0 − Re )+ k1 ( R0 − Re )

3

] t$ − [c R& + c R& ] t$ E

z

0

e

1

3 e

E

z

(7.3.2)

where R0 is the free length of the spring (unloaded wheel radius).

Tractive Forces The Magic Formula Tyre Model requires as inputs the longitudinal slip, the lateral slip and the normal force on the tire. The normal force is easily obtained from the vertical force model as discussed previously. The longitudinal and lateral slip angles are somewhat more difficult to obtain and are discussed below. The output of the Magic Tyre Formula consists of a system of forces and moments at the center of the contact patch. In order to keep the tire model to vehicle model interface as general as possible it is desirable to eliminate any references to point(s) of contact. It is preferable to return a system of forces and moments applied at the wheel center. This also eliminates the need to determine the virtual displacements associated with any points of contact. It is trivial to relocate a system of forces and moments; the forces and moments are modified as follows: E

Fwc = E Fcp

121

(7.3.3)

E

M wc = E M cp + E rCP/ WC ×E Fcp

(7.3.4)

where a WC subscript indicates the force or moment as applied at the wheel center and a CP subscript indicates the force or moment as applied at the contact point. Lateral Slip and Longitudinal Slip The SAE definition of longitudinal slip velocity is ω-ω0 where ω is the actual angular velocity of the tire and ω0 is the angular velocity of a free-rolling tire moving with the same linear velocity as the driven or braked tire. The longitudinal slip

percentage is

defined as the ratio of the slip velocity to the free rolling angular velocity:

κ=

ω −1 ω0

(7.3.5)

The SAE definition works well for forward velocities but breaks down when negative velocities are considered. The desired result is for the longitudinal force to have the same sign as the tractive force under all conditions of braking and acceleration for both forward motion and backward motion of the vehicle. The SAE definition can be modified to product the proper results by taking the magnitude of the angular velocity in the denominator:

κ=

ω −1 ω0

(7.3.6)

or written in terms of linear velocities and extending to negative velocities,

κ= −

Vx − Vr V = − sx Vx Vx

(7.3.7)

where Vsx = Re (ω 0 − ω ) is also referred to as the longitudinal slip velocity, Vr = Re ω is the forward speed of rolling, V x = Re ω 0 is the actual forward velocity of the tire and Re is

122

the effective radius of the tire. The relationships between the Vi are

Vr

Vsy V

Vs

shown in Figure 7.2. Table 7.1

α

shows that the preceding expression produces the desired sign for κ Vsx

Vx

under all conditions of forward and Figure 7.2: The Relationship between the Components of the Tire Velocity Vector

reverse motion. A problem arises due to the

singularity in the slip percentage that occurs at Vx=ω0=0. This singularity is encountered when the linear velocity of the tire goes to zero. The solution to this problem is to derive a first order differential equation for longitudinal slip following the example of Bernard and Clover1 [Bernard, 1995],

κ& +

Vx V κ = − sx B B

(7.3.8)

where Vx is the longitudinal component of the wheel velocity, and B is an experimentally determined parameter with units of distance which characterizes the first order lag. Note that the steady state solution to the preceding equation is identical to the SAE definition of slip. Note that the differential equation has good behavior for small values of Vx including Vx=0.

1

Bernard and Clover’s definition of longitudinal slip is the negative of the SAE definition used by Pacejka’s Magic Formula Tyre Model. Bernard and Clover’s results have been modified to use the SAE definition of slip.

123

Table 7.1 - Desired Longitudinal Force Sign and Sign of Longitudinal Slip

Vx > 0 (Forward Motion)

Vx < 0 (Reverse Motion)

Vr < Vx

Vr > Vx

(Braking) Desired: Fx < 0 Sign of κ: Desired: Fx > 0 Sign of κ: +

(Acceleration) Desired: Fx > 0 Sign of κ: + Desired: Fx < 0 Sign of κ: -

Bernard and Clover found that oscillations occur in the lateral and longitudinal slip at low velocities. To eliminate these oscillations it is necessary to include a damping term which is only activated when the velocity is below a certain threshold (0.15 m/s is suggested) and changes sign from one integration time step to the next. The form of the damping used by Bernard and Clover is

Fdamping ,x = 2ζ Vx

Fz Cs gB

(7.3.9)

The term Fz/g can be viewed as the portion of the vehicle mass supported on the wheel in question, Cs is the longitudinal stiffness, ζ is the damping coefficient and Vx is the longitudinal velocity of the tire. The lateral slip angle can be treated with a similar approach. The SAE definition for slip angle is the angle between the vector defined by the intersection of the wheel plane and the road plane and the velocity vector of the center of the contact patch. It is necessary to modify the SAE definition to handle negative velocities as was done for the longitudinal slip. With this modification the slip angle is given by tan α = −

Vsy Vx

124

(7.3.10)

where u and v are the lateral and longitudinal components of the wheel velocity vector. Modifying Bernard and Clover’s derivation for slip angle yields

V V d (tan α )+ x tan α = − sy dt b b

(7.3.11)

Note that the steady state solution matches the modified SAE definition of slip angle. The parameter b is the relaxation length for the tire which controls slip angle lag. A damping force is also applied to the wheel in the lateral direction to eliminate oscillations in the lateral velocity which occur at small velocities.

Fdamping , y = 2ζ Vsy

Fz Cα gb

(7.3.12)

Magic Formula The Magic Formula Tyre Model 2 is a complete tire model in that it allows determination of all six of the forces and moments generated by the tire: lateral force, longitudinal force, normal force, rolling resistance, overturning moment and self-aligning torque. The normal force model specified by the Magic Formula Tyre model consists of a simple linear spring and linear damper. This component of the model has been modified to include nonlinearities as discussed in the preceding section. The driving/braking moment is modeled as part of the power train and braking submodel. The equations for each of the forces or moments are based on one of the two Magic Formulas:

2

The version of the Magic Formula Tyre model used here is the static version of the Delft Tyre ’97 model which is designed to handle the combined slip case. See Pacejka (1997) for a complete description of the model. The previous version of the Magic Formula Tyre model is described in Pacejka (1992) and in Genta’s text. Additional information on even earlier forms of the model can be found in Bakker’s papers.

125

Y = y + SV

[

{

}]

(7.3.13)

{

}]

(7.3.14)

y = D sin C arctan Bx − E( Bx − arctan ( Bx )) x = X + SH Y = y + SV

[

y = D cos C arctan Bx − E( Bx − arctan ( Bx )) x = X + SH

In both of the formulas the parameter D is the maximum value the force or moment apart from the small effect due to the Sv term. For the sin() form of the formula the product BCD gives the slope of the curve at σ + S h = 0 which, in the case of the lateral force model, corresponds to the initial cornering stiffness. Sv and Sh are introduced to allow for non-zero forces and moments at zero slip. This can occur due to asymmetries in the tire’s construction (ply steer, conicity, etc.). The parameter C is called the shape factor and limits the range of the arguments in the sin() function. This leaves the factor B to control the slope of the curve at the origin and hence it is referred to as the stiffness factor. For the cos() form of the formula the product BC controls the breadth of the peak. The parameter C controls the value of the asymptote of the function as x goes to positive or negative infinity and thus acts to shape the flanks of the curve. This leaves the parameter B to control the peak width. As the Magic Tyre Formula has evolved it has become necessary to express the parameters B through E as functions of other variables such as normal load ( Fz) and camber angle ( γ ) in order to obtain the required levels of accuracy (see Pacejka, 1997 for details). The combined slip case is handled by an additional set of formulas, also based on the form of the ‘magic’ formula, which take as

126

arguments the forces and moments for the pure slip cases. Equations are also given for rolling resistance force, overturning moment and normal load determination.

7.4 Generalized Forces and Moments Once the forces and moments applied at the wheel center by the tire have been determined the generalized forces and the generalized moments can be determined. The calculation of the generalized forces is slightly different for the front and for the rear tires due to the differences in the suspension geometry.

Front Tires The position of the center of the wheel for the front tires is specified relative to the spindle/steering knuckle assemble. The position of the wheel center is E

rWC/ E = E rSP/ E + [E CSP ]SP rWC/SP

(7.4.15)

where the SP subscript indicates the spindle and WC indicates the wheel center. Computation of the generalized forces gives the following results: δErWC/ E = (1 0 0)δxSP + (0 1 0)δySP + (0 0 1)δzSP +

 ∂[E CSP ]SP  rWC/SP δβSP,i ∂βSP,i 

∑  i

(7.4.16)

T δW = E Fspring δErWC / E

(7.4.17)

T QxSP = E Fspring (1 0 0)

(7.4.18)

T Q ySP = E Fspring ( 0 1 0)

(7.4.19)

T QzSP = E Fspring ( 0 1 0)

(7.4.20)

and

127

 ∂[E CSP ]SP T   rWC /SP QβSP,i = E Fspring  ∂βSP,i 

(7.4.21)

The virtual work performed by an applied moment can be shown 3 to be of the form δW =2 E M TWC [E G SP ]δβ

(7.4.22)

The component of the moment generated by the tire which lies along the axis of rotation of the wheel is responsible for the wheel’s rotation. The remaining components are transferred to the spindle via the wheel bearings. Dividing the moment into these components gives E

M ROT = (0 1 0)[E CSP ]E M WC

(7.4.23)

M SP = E M WC − E M ROT

(7.4.24)

and E

The generalized forces are thus T QβSP,0 =2 E M SP [E G SP ](1 0 0 0)

(7.4.25)

T QβSP,1 =2 E M SP [E G SP ](0 1 0 0)

(7.4.26)

T QβSP,2 =2 E M SP [E G SP ](0 0 1 0)

(7.4.27)

T QβSP,3 =2 E M SP [E G SP ](0 0 0 1)

(7.4.28)

QΩ = E M TROT

(7.4.29)

T

T

T

T

where

3

See Nikravesh, p.290 for a similar problem. Extending Nikravesh’s derivation to a force couple gives the desired result.

128

 − βSP,1 [E G SP ]= − βSP,2   − βSP,3

βSP,0

− βSP,3

βSP,3 − βSP,2

βSP,0 βSP,1

βSP,2   − βSP,1  βSP,0  

(7.4.30)

Rear Tires The position of the center of the wheel for the rear tires is specified relative to the rear suspension coordinate system RS which is in turn located relative to the sprung mass. The position of the wheel center is E

rWC / E = E rSM / E +

[C ]

SM

E

SM

rpiv /SM +

[ C ]( E

RS

RS

rRS/ piv +

RS

rWC/ RS

)

(7.4.31)

where WC indicates the wheel center. Computation of the generalized forces gives δErWC/ E = (1 0 0)δxSM + (0 1 0)δySM + (0 0 1)δzSM +

 ∂[E CSM ]SM  rpiv /SM δβSM,i + ∂βSM, i 

∑  i

δW =

(

E

 ∂[E C RS ]  ∂βRS,i 

(

∑  i

RS

rpiv / RS +

)

Fspring + E Fdamper ⋅δE rWC / E

RS

)

rWC/ RS δβRS,i

(7.4.32)

(7.4.33)

The generalized forces are simply the coefficients of the virtual displacements in the expression for the virtual work: T QxSM = E Fspring (1 0 0)

(7.4.34)

T Q ySM = E Fspring ( 0 1 0)

(7.4.35)

T QzSM = E Fspring (0 0 1)

(7.4.36)

∂[E CSM ] SM T QβSM,i = E Fspring   rpiv /SM  ∂βSM,i 

(7.4.37)

 ∂[E C RS ] T QβRS,i = E Fspring    ∂βRS,i 

(

RS

129

rRS/ piv +

RS

rWC/ RS

)

(7.4.38)

Dividing the moment into components gives E

M ROT = (0 1 0)[E C RS ]E M WC

(7.4.39)

M SP = E M WC − E M ROT

(7.4.40)

E

Computing the generalized forces gives

QβRS,0 =2 E M TWC [E G RS ](1 0 0 0)

(7.4.41)

QβRS,1 =2 E M TWC [E G RS ](0 1 0 0)

(7.4.42)

QβRS,2 =2 E M TWC [E G RS ](0 0 1 0)

(7.4.43)

QβRS,3 =2 E M TWC [E G RS ](0 0 0 1)

(7.4.44)

QΩ = E M TROT

(7.4.45)

T

T

T

T

where

 − βRS,1 [E G RS ]= − βRS,2   − βRS,3

βRS,0

− βRS,3

βRS,3 − βRS,2

βRS,0 βRS,1

130

βRS,2   − βRS,1  βRS,0  

(7.4.46)

8 Equations of Motion - Driver Model 8.1 Introduction The driver model is a critical component of the simulation. It is responsible for controlling the steering, acceleration and braking inputs to the vehicle. In a racing situation, lap time is strongly dependent on the ability of the driver to operate the vehicle at or near the limits of handling and performance. This skill must be reflected in the driver model. Therefore, the degree of skill with which the computer performs the driving task plays a significant role in the accuracy and usefulness of the simulation. There are several classes of driver models which appear frequently in the literature, either alone, or in combination with one another. These include quasi-linear compensatory controllers, pursuit mode controllers and precognitive controllers. In general, quasi-linear compensatory controllers, which utilize feedback loops based on the current vehicle position and orientation, are unable to perform adequate lane keeping control due to the lag time in the vehicle response. Pursuit controllers attempt to model the driver’s ability to see the roadway ahead and can thus initiate control input prior to reaching a turn or arriving at a braking point. Precognitive controllers model the human driver’s ability to execute learned maneuvers on command. This is typically implemented by specifying a steering wheel angle directly without regard to feedback quantities. A number of these

131

driver control approaches were tried before the final design was chosen for both the steering control case and for the throttle and brake control case.

8.2 Steering Control The first several steering control attempts utilized pursuit type controllers. The most successful of these was based on the optimal preview-follower theory of Guo and Guan (Guo, 1993). This controller worked reasonably well under steady state velocity conditions, low lateral accelerations and on flat road surfaces. The weakness of this controller, however, is it’s dependence on knowing the transfer function which relates the lateral acceleration of the vehicle to the steering angle. While this transfer functions can be readily determined for a simple linear vehicle model with ideal tires, it must be generated computationally for more complex vehicle models containing significant nonlinearities. The lateral response transfer function for both the linear model and the nonlinear model is strongly dependent on vehicle speed, lateral acceleration and road bank angle. It was deemed impractical to characterize the lateral response transfer function for the nonlinear vehicle model over the broad range of operating conditions necessary for the simulation. Characterization of the vehicle is even more difficult during the optimal design process because changes made to the vehicle suspension setup during the optimization process can produce significant changes in the response functions. This necessitates a re-tuning of the driver controller during the optimization process. The inherent problems with compensatory and pursuit controllers discussed above made it necessary to try a different approach to the driver control. The resulting controller

132

Figure 8.1 - The Driver Path for the Kenley, NC Race Track

fits into the precognitive controller category in that the steering angle is specified as a function of the vehicle’s position along the roadway and is specified without regard to feedback quantities. On the other hand, the steering angle vs. road position curve is obtained via minimization of a cost function based on the accumulated lateral position error. In this sense the controller can be thought of as an optimal controller which can adapt to changes in the vehicle setup and which can operate at the limits of adhesion. The details of the driver path definition, of the cost function computation, and of the optimization process are discussed below.

Driver Path Definition The driver path is specified relative to the centerline of the road. The lateral offset of the path is specified at the beginning, middle and end of each road segment. A quadratic

133

polynomial is fitted through the three points to form a smooth path through the segment. The current version of the driver path formulation does not attempt to match tangents at the road segment boundaries. The driver path segments are adjusted by eye to obtain reasonable continuity of the driver path tangent vector between road segments. The driver path, indicated by the dark line, for the Kenley, NC race track is shown in Figure 8.1. The positioning of the driver path is based on the experience of the user of the software and should approximate the line followed by real drivers when racing on the real track. At present no optimization of the driver path is performed by the computer. The steering curve is described as a series of points which consist of the steering angle and the position of the point along the track. The position coordinate is specified as a distance from the start of the track as measured along the centerline of the road. The lower bound on position is zero and the upper bound on the position is the length of the track as measured along the centerline. The steering angle between points is determined via linear interpolation. A sample steering profile utilizing ten data points is shown in Figure 8.2.

Steering Profile Optimization and Cost Function Computation As mentioned before, an optimal steering profile is obtained via the minimization of a cost function based on lateral position error. It is necessary to include the parameters which describe the steering profile in the optimization process, along with those parameters of the vehicle model which are being optimized, since the handling characteristics of the vehicle may change as part of the optimization. This has the unfortunate side effect of significantly increasing the dimension of the optimal design 134

45.0 40.0 35.0

Angle (degrees)

30.0 25.0 20.0 15.0 10.0 5.0 0.0 -5.0 -10.0 0

200

400

600

800

1000

1200

1400

1600

1800

2000

Track Position (ft) Steering Wheel Angle

Unconstrained SW Profile Point

Constrained SW Profile Point

Figure 8.2 - The Steering Profile for the Kenley, NC Race Track

search space and of slowing the optimization process. On the other hand, it allows the simulation to consistently drive the vehicle to it’s limits. Ordinarily, two parameters are used to describe each steering profile point: position and steering wheel angle. For the ten point steering profile shown in Figure 8.2, this leads to a 20 parameter optimization space. The size of the optimization space is further increased by the addition of the speed profile parameters (discussed below) and by the addition of vehicle design parameters. The number of parameters being optimized can rapidly become unwieldy for any but the simplest of road courses. Fortunately, the number of parameters can be significantly reduced by taking advantage of symmetry and/or

135

periodicity in the steering profile and also by fixing the position of those steering points for which optimization of both the position and the steering wheel angle is unnecessary. The cost function used to optimize the steering profile is computed by integrating the magnitude of the instantaneous error in lateral position with respect to position along the track and then dividing by the length of the track.

C.F.=

1 L y veh ( x ) − y path ( x ) dx L ∫0

(8.2.1)

The resulting value is the average lateral position error of the vehicle measured over the entire lap. The lateral position of the vehicle is determined by projecting the vehicle’s reference point (usually taken to be a point near the geometric center of the vehicle) onto the road surface and then finding the distance between the projected point and the road centerline along a line perpendicular to the road centerline’s tangent vector. The lateral position of the driver path is determined by evaluating the interpolating polynomial which defines the driver path at the position corresponding to the point on the road centerline which was just determined. The computation of the cost function can be facilitated by integrating it along with the equations of motion. Integrating the cost function in the time domain requires that the variable of integration in Equation 8.2.1 be changed: dx = v veh dt C.F.=

⇒

dx = v veh dt

1 T y veh ( t ) − y path ( t ) v veh ( t )dt L ∫0

136

(8.2.2)

(8.2.3)

where vveh is the velocity of the vehicle at the time (or position) in question. Note that integrating the path error in time without the velocity weighting term would lead to a cost function which de-emphasizes errors in the sections of the track which are traversed quickly. An even weighting is far more desirable. In the event that the vehicle crashes before the completing the desired number of laps or before reaching the end of the course, a penalty is computed and added to the cost function value as computed up to the time of the crash. In designing the function for the penalty term computation there are two desirable characteristics which should be satisfied if possible. The first characteristic is that the penalty for a crash near the end of the road course should carry less weight than the penalty for a crash near the beginning of the course. This property aids the optimization of the steering profile by encouraging the optimizer to move towards a path which results in the vehicle completing the entire road course. The second desirable characteristic of the penalty term is that, in the worst case, it should be at least as large as the maximum possible cost function value which could be obtained by a vehicle which successfully complete the road course. A penalty function which satisfies both of these conditions can be obtained by integrating the maximum possible lateral position error for the remainder of the distance around the track. The penalty term is computed by integrating the function

P.F.=

xend

∫

xcurr

ymax ( x ) dx

137

(8.2.4)

where xcurr is the current position (the position at the time of the crash), xend is position at the end of the road course and ymax is the maximum distance between the driver path and the edge of the road as defined below:

 1 w + y path ( x ), y max ( x ) =  21  2 w − y path ( x ),

y path ( x ) ≥ 0

y path ( x ) < 0

(8.2.5)

8.3 Speed Control The speed control system is similar to the steering control system in that a prescribed velocity profile is defined which specifies the desired vehicle velocity as a function of position along the track. Unlike the steering control system the profile is not used as a direct input to the vehicle model. Instead, the goal of the speed controller is to modulate the throttle and brake to follow the speed profile as precisely as possible. The controller utilizes a single point preview strategy to anticipate changes in the vehicle velocity. The use of preview allows the controller to begin responding to a rapid change in the velocity profile before the vehicle actually arrives at the point where the change occurs, thus reducing error in tracking the velocity profile. The controller also incorporates traction control and anti-lock brake control features which operate by limiting the longitudinal slip at the tires. This is done in order to prevent wheel lockup or excessive wheel spin which, if left uncontrolled, can cause numerical problems with the integration process. The speed profile is specified in the same manner as the steering profile: A series of speed/position data points defines the basic curve and linear interpolation is used to obtain

138

s Vdes

P( s)

Veff

Gerr

Aeff

Ain

C( s)

Flong

V ( s)

Aveh 1

Vveh s

Figure 8.3 - The Driver Speed Controller Block Diagram

speed values between the data points. The parameters describing the speed profile are typically included in the optimization process, the goal being to minimize the lap time by maximizing the vehicle speed. Note that, since the velocity profile can change during the optimization process, it may be necessary to dynamically update the initial conditions used for each simulation. The block diagram for the speed control algorithm is shown in Figure 8.3. The form of the control was generated by considering the types of information available to the speed controller and the types of input (external inputs and feedbacks) which could be applied to the vehicle model. The available inputs and feedbacks were then combined to produce the desired vehicle model control inputs. The detailed discussion of the driver control algorithm, presented below, follows the logical progression used in its development: identifying the available inputs and feedback quantities, identifying and generating the control input to the vehicle model, and finally, the adding the preview functionality.

139

The sole external input consists of the prescribed velocity profile (designed by Vdes in the block diagram). By differentiating the velocity profile with respect to time the desired acceleration at any point along the track may also be determined. Since the velocity profile is specified as a function of track position it is necessary to apply the chain rule for differentiation to find the acceleration at a particular point along the track. d Vdes ( x ) dt dx d = Vdes ( x ) dt dx  d = Vdes ( x ) Vdes ( x )   dx

Ades ( x ) =

(8.3.6)

The desired acceleration is used as the primary input to the velocity controller. Under ideal circumstances it would be the only necessary input since tracking the desired acceleration exactly should yield the desired velocity curve. In practice it is necessary to included an additional feedback loop to provide some error correction ability. There are a number of potentially useful feedback terms generated by the vehicle model including vehicle speed, wheels speed, longitudinal slip, lateral slip and so on. The most useful of these is, of course, the longitudinal velocity of the vehicle, designated by Vveh in the block diagram. The vehicle’s longitudinal velocity is subtracted from the previewed desired longitudinal velocity to generate an error value which is multiplied by a gain Gerr (which has units of s-1). The resulting value is then added to the previewed acceleration value (the preview block is discussed in greater detail below) and provided as a command input to the driver dynamics control block.

140

Driver Dynamics Block The dynamics of the driver are represented by the control block C(s). This block converts the command acceleration to a force to be applied to the vehicle via the tires and implements the traction control and anti-lock brake features. Unlike many of the driver control models found in the literature, there is no need to model the lags and the frequency response limitations of a human driver in this application. In fact, modeling these delays would most likely degrade the tracking performance of the controller. In light of these considerations, the conversion portion of the driver control block consists of a simple gain, equal to the effective mass of the vehicle, which converts the acceleration command input to a force to be applied to the vehicle via the rear tires, in the case of acceleration, or via all four tires, in the case of braking. Facc _ brk ( t ) = M eff ⋅ Ain ( t )

(8.3.7)

The manner in which the force is applied to the vehicle depends on whether the force represents an accelerating force or whether it represents a braking force. In the acceleration case, the longitudinal force is divided equally between the rear wheels (since the rear differential is locked in a Legends car) and converted to an applied moment by dividing by the current wheel radius. In the braking case, the force is apportioned between the front and rear of the car using a fixed brake bias constant and then divided equally between the left and right wheels. Again, the force is converted to an applied moment by dividing by the wheel radius. The traction control portion of the driver control block is somewhat more complex. The goal of the traction control algorithm is to limit wheel spin under conditions

141

of excessive drive torque. This is accomplished by reducing the drive torque in proportion to the degree of excessive wheel slip. To design a useful controller it is first necessary to understand the physical significance of the longitudinal slip and to determine the possible range of longitudinal slip values. The range of the longitudinal slip variable under conditions of acceleration can be determined by considering the steady state solution to the differential equation for longitudinal slip introduced in the chapter on tire modeling:

κ& +

Vx V κ = − sx B B

(8.3.8)

which has the steady state solution

κ=−

Vx − Vr Vx

(8.3.9)

where Vx is the actual longitudinal velocity of the wheel center and Vr is the velocity of the wheel center if it were in a free rolling state (i.e. Vr = Rω where R and ω are the current wheel radius and angular velocity). Under conditions of acceleration (either forward or reverse) the magnitude of the rolling velocity Vr is greater than the magnitude of the actual velocity Vx. This condition produces values of κ greater than zero. If the angular velocity of the wheel is zero (locked wheel under braking) κ is negative one. Under conditions of extreme acceleration where Vr is much greater than Vx, κ approaches positive infinity. A value of κ equal to one is obtained when Vr = 2Vx or the wheel is spinning at twice the free rolling angular velocity.

142

1.1

Acceleration Gain Factor (unitless)

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Longitudinal Slip (unitless) Gtc = 0.0 Gtc = 1.0

Gtc = 0.05 Gtc = 1.5

Gtc = 0.1 Gtc = 2.0

Gtc = 0.25

Gtc = 0.5

Figure 8.4 - The Effect of the Traction Control Gain Parameter on the Acceleration

To control wheel spin under acceleration it is desirable to gradually reduce the magnitude of the acceleration command input, and thus the applied torque, as the amount of wheel spin increases. The reduction should begin at the desired slip threshold and increase until the applied torque goes to zero at infinite slip. A gain factor which can be applied to the command acceleration, thus reducing the drive torque, and which has the desired characteristics is

Gaccel

10 . ,κ + Coffset ≤ 10 .   1 = ,κ + Coffset > 10 .  1 + GTC (κ + Coffset − 1) 

(8.3.10)

The constants GTC and Coffset control the rate of decrease of the command acceleration and the threshold at which the traction control begins to take effect respectively. Figure 8.4 shows a plot of the Gaccel as a function of longitudinal slip for several values of GTC.

143

Vehicle Dynamics Block Ideally, the tires would immediately generate the desired longitudinal force upon application of an appropriate driving torque. In reality this doesn’t happen for two reasons. The first reason is that the longitudinal force generated by a tire lags the application of a driving or braking torque. This lag is modeled by the first order longitudinal slip equation in the tire model and the lag is characterized by the longitudinal relaxation length for the tire. These lags are represented in the block diagram of the speed controller by the vehicle dynamics block V(s):

V (s) =

1 τ lag s + 1

(8.3.11)

where τ lag is the time constant associated with slip equations. Note that the time lag may vary with vehicle velocity. The second reason that the longitudinal force may not track the applied torque exactly is that the applied torque may exceed the tractive capabilities of the tire. The traction control algorithm does nothing to alleviate this problem since all it does it to limit wheel spin under conditions of excessive acceleration and it does not actually enhance the tractive capabilities of the tires. These two effects are not modeled in the controller block diagram.

Preview Compensation Block The preview block is used to partially compensate for the lags inherent in the vehicle response discussed above. As indicated by the name, the function of the preview block is to generate an input to the remainder of the vehicle speed controller which is based on profile information from the section of road ahead of the vehicle. This enables

144

the controller to anticipate changes in the vehicle’s velocity and thus reduce lags in the response. There are number of ways to model the driver’s ability to see the road (or in this case the velocity profile) ahead of the vehicle. The majority of the models utilize a weighted average of the road profile or of the velocity profile information directly ahead of the vehicle. The weighted average usually takes a form similar to τ2

Veff ( t ) =

∫ w(τ )V ( t + τ )dτ ∫ w(τ )dτ τ1

des

(8.3.12)

τ2

τ1

where Vdes(t) is the prescribed velocity profile, Veff(t) is the effective velocity input provided to the control algorithm (the previewed velocity) and w(τ) is a weighting function. The integration limits define the boundaries of the preview interval. The simplest type of preview function to implement is the so-called single point preview. Single point preview is used in this application. The single point preview form is obtained by inserting the weighting function

(

w ( τ ) = δ t + Tp

)

(8.3.13)

where δis the Dirac delta function, into the preceding equation. The result is

(

Veff ( t ) = Vdes t + Tp

)

(8.3.14)

where Tp is the preview time. The transfer function for the preview block can be obtained by computing the Laplace transform of the preceding result.

{ (

Veff ( s) = L Vdes t + Tp

)}= e L{V ( t )}= e Tp s

145

Vdes ( s)

Tp s

des

(8.3.15)

V des

V eff

P(s)

F (s)

V veh

Figure 8.5 - The Simplified Preview-Follower Control System

The resulting transfer function can be expanded in a Taylor’s series to get a polynomial form which is used in later computations.

P( s) =

Veff ( s) Vdes ( t )

=e

Tp s

= 1 + Tp s +

Tp2 2

s2 + K

(8.3.16)

= 1 + P1s + P2 s 2 + K At this point, the form of the preview function has been chosen. It is still necessary to tune the preview block to provide the best velocity tracking performance by determining the proper value for Tp. In order to obtain an optimal controller, it is necessary to consider the dynamics of the combined controller-vehicle system. The dynamics of the system can be depicted using a simplified block diagram like the one shown in Figure 8.5. The dynamic response of the preview function is represented by the P(s) (preview) block and the dynamic response of the combined vehicle and control network is represented by the F(s) (follower) block. For an ideal control system the product of the preview transfer function and of the follower transfer function should be unity for all values of s. This would produce an output equal to the input (i.e. perfect tracking of the velocity profile). In reality, it is sufficient for this equality to be satisfied within a particular low frequency range (usually below 5-10 Hz) [Guo, 1993]. P( s) ⋅ F ( s) ≈1

146

(8.3.17)

If the inverse of the follower transfer function, which represents the combined controller response and vehicle response, is given by 2 3 F − 1( s) = 1 + Ps 1 + P2 s + P3s + K

(8.3.18)

then the product of the preview function and the follower function (which represents the complete system) is unity within the low frequency range and perfect control is achieved. In practice it has been shown that matching the coefficients beyond third order can lead to controller stability problems [Guo, 1993]. On the other hand, matching an insufficient number of coefficients can lead to significant tracking errors. The transfer function for the controller shown in Figure 8.3 can be shown to be

 1   ( Gerr + s)V ( s) vveh ( s)  s = F ( s) = veff ( s)  1 1 +   Gerr V ( s)  s

(8.3.19)

where V(s) is the block representing the vehicle dynamics and where the effects of traction control and anti-lock brake control on the system are ignored. Setting the product of the preview transfer function and the follower transfer function equal to one gives the following result

 1  1   ( Gerr + s)V ( s)P( s) = 1 +   Gerr V ( s)  s  s

(8.3.20)

Inserting the Taylor’s series expansion of the preview transfer function (keeping only constants and terms linear in s) derived earlier and canceling terms gives

[

]

V ( s) 1 + Gerr Tp + Tp s = 1

147

(8.3.21)

The transfer function used to model the vehicle response to acceleration commands was assumed to be a simple first order lag:

V ( s)=

1 1 + τ lag s

(8.3.22)

Inserting this result into the preceding equation and collecting terms based on the order of s gives two equations: Tp = τ lag

(8.3.23)

Gerr Tp = 0

(8.3.24)

The first equation provides a means of selecting the best preview time for a given vehicle once the vehicle’s response has been characterized. The second equation indicates that Gerr should be zero (since Tp is non-zero according to the first equation). This effectively removes the velocity feedback loop from the controller. While this might work for a vehicle whose response matches the vehicle dynamics model V(s) exactly, the real vehicle requires velocity error feedback to correct deviations from the desired response. The deviations are primarily due to limitations of the tires at high slip values and to the action of the traction control and anti-lock braking control algorithms at high acceleration values. Even though Gerr can’t be set to zero it is desirable to satisfy the condition as nearly as possible by keeping the product of the preview time and the gain reasonably small (so that the acceleration command generated by the feedback loop is small compared to the acceleration command from the prescribed velocity profile). The preview time can be determined once an estimate for the lag constant τ lag has been computed. The lag constant can be estimated given the relaxation length of the tire

148

and the velocity of the vehicle. The differential equation for the longitudinal slip of the tire is of the form

Vx V κ = − sx B B

κ& +

(8.3.25)

where the time constant for the exponential solution is

τ lag =

B Vx

(8.3.26)

For a typical longitudinal velocity on the order of 100 ft/sec and a typical relaxation length of 0.3 feet the lag time is found to be τ lag = 0.003 sec

(8.3.27)

Inspection of simulation results seems to back up this result with lag time values on the order of 0.1 seconds or less. The gain Gerr controls the responsiveness of the controller to errors in the vehicle velocity. Selecting a large gain produces rapid responses but can interfere with the desired acceleration feed forward. Selecting a gain which is too small causes long delays in correcting velocity errors. Experimentation has shown that taking Gerr ≈ 0.5 works well when simulating the NCSU Legends car. Note that the product of the of the velocity error gain and the preview time is small (on the order of 0.01) which comes close to satisfying condition 8.3.24.

149

9 Results, Conclusions and Recommendations 9.1 Introduction The ultimate goal of the vehicle dynamics research program at NCSU is to demonstrate the applicability of optimal design techniques used in other automotive applications (e.g. valvetrain and cam design [Etheridge, 1998 and Kim, 1990]) to vehicle dynamics. While complete success has not yet been achieved, a great deal of progress has been made. A vehicle model with an intermediate number of DOF has been developed and a computer program has been written to solve the equations of motion. The computer simulation has demonstrated the necessary computational efficiency to allow optimal design to be performed. The vehicle model and the associated computer simulation have been partially validated by modeling the NCSU Legends race car. The following sections discuss the measurement of the NCSU Legends Car, present the vehicle data used in the simulation, discuss the chassis setup procedure applied to the model, and finally, discuss the simulation results.

9.2 Measurement Process and Model Data The first step of the simulation and/or optimization process is to determine the geometric parameters and physical parameters which describe the system being modeled. The system in this case consists of the vehicle itself, the tires, the road surface and the

150

driver. The measurement processes used to obtain the data and the data itself are presented in the following sections for each of the components of the system.

Vehicle Data The vehicle data primarily consists of the physical data describing the positions of the various joints connecting the components which make up the vehicle and the mass and inertia properties of those components. The vehicle measurement process is described below and is followed by the vehicle data. The first step in the vehicle measurement process is to locate the origins of the various body fixed coordinate systems. The use of centroidal body fixed coordinate systems in the derivation of the equations of motion requires that the origin of the each of the coordinate systems used to measure the car be located at the center of gravity of the respective bodies. Thus, it is necessary to compute the location of the center of gravity of each model component prior to measuring the vehicle. The orientation of the body fixed coordinate systems are subject to a few restrictions which result from additional conditions imposed on the system during the

Table 9.1 - Front Suspension Geometric Data

Parameter Description

Left Front (ft)

Right Front (ft)

Wheel Center Location Lower Ball Joint Location Upper Ball Joint Location Steering Knuckle Location Control Arm Spring Mounts Control Arm Damper Mounts Upper Control Arm Length Lower Control Arm Length

{0, 0.0833, 0 } {0, -0.375, -0.315 } {-0.015, -0.41, 0.3 } {0.325, -0.323, -0.305} {0.125, 1.15, 0 } {0.125, 1.15, 0 } 0.677 1.33

{0, -0.0833, 0 } {0, 0.375, -0.406 } {-0.015, 0.41, 0.3 } {0.325, 0.281, -0.333} {-0.125, 1.13, 0 } {-0.125, 1.13, 0 } 0.688 1.35

151

derivation of the equations of motion. The x-axis of the sprung mass coordinate system should be parallel to the longitudinal plane of symmetry of the vehicle. This requirement provides the lateral symmetry required by the steering system model. The x-axis need not be coincident with that plane. The y-axis of the unsprung mass coordinate systems should be parallel to the axes of rotation of the wheels. This requirement is necessary to preserve the rotational symmetry of the inertia tensor of the wheels in the unsprung mass body fixed coordinate systems. Based on these considerations, the coordinate systems used for the four masses are defined as follows. The sprung mass x-axis points forward from the center of gravity parallel to the plane of lateral symmetry. The y-axis points out the driver’s side door and the z-axis points up. The remaining three unsprung mass coordinate systems are oriented

Table 9.2 - Rear Suspension Geometric Data, Spring Data and Damper Data

Parameter Description

Value (ft)

Suspension Link Lengths: Left Trailing Link Center Trailing Link Right Trailing Link Panhard Rod Suspension Link Mounting Locations: Left Trailing Link Center Trailing Link Right Trailing Link Panhard Rod Spring and Damper Mounting Locations: Left Side Spring Left Side Damper Right Side Spring Right Side Damper Wheel Center Locations: Left Side Right Side

1.23 0.775 1.2 2.19 { -0.0167, 1.29, -0.271 } { -0.133, -0.0933, 0.354 } { 0.00833, -0.883, -0.312 } { -0.192, 1.23, -0.283 } { 0.292, 1.67, -0.275 } { 0.292, 1.67, -0.275 } { 0.292, -1.26, -0.282 } { 0.292, -1.26, -0.282 } { 0, 2.33, 0 } { 0, -1.83, 0 }

152

in approximately the same manner. The y-axis for each of the unsprung mass coordinate systems is parallel to the axis of rotation of the wheels. With the location of the origin and the orientation of the coordinate axes having been determined for each of the four body fixed coordinate systems, the measurement process can be begun. Since it is rarely convenient to measure the locations of points on

Table 9.3 - Sprung Mass Geometric Data

Parameter Description

Value (ft)

Spring Mounting Locations: Left Front (Upper) Right Front (Upper) Left Rear (Upper) Right Rear (Upper) Damper Mounting Locations: Left Front (Upper) Right Front (Upper) Left Rear (Upper) Right Rear (Upper) Control Arm Coordinate System Origin: Left Front (Upper) Left Front (Lower) Right Front (Upper) Right Front (Lower) Control Arm Axes of Rotation: Left Front (Upper) Left Front (Lower) Right Front (Upper) Right Front (Lower) Rear Suspension Mounting Points: Left Trailing Link Center Trailing Link Right Trailing Link Panhard Rod Reference Point Locations: Left Front Passenger Compartment Frame Corner Right Front Passenger Compartment Frame Corner Left Rear Passenger Compartment Frame Corner Right Rear Passenger Compartment Frame Corner Vehicle Center

153

{ 3.31, 1.18, 0.771 } { 3.29, -1.29, 0.792 } { -3.32, 1.08, 1 } { -3.29, -1.19, 1 } { 3.31, 1.18, 0.771 } { 3.29, -1.29, 0.792 } { -3.32, 1.08, 1 } { -3.29, -1.19, 1 } { 3.13, 0.9, 0.344 } { 3.1, 0.26, -0.281 } { 3.13, -1.04, 0.354 } { 3.1, -0.385, -0.271 } { 1, 0, 0 } { 1, 0, 0 } { -1, 0, 0 } { -1, 0, 0 } { -2.41, 1.25, -0.167 } { -2.95, -0.354, 0.438 } { -2.39, -1.29, -0.167 } { -4.05, -1.22, 0.0417 } { 1.73, 0.750, -0.4583 } { 1.73, -0.883, -0.4583 } { -2.73, 1.18, -0.4583 } { -2.73, -1.31, -0.4583 } { 0, 0, 0 }

Table 9.4 - NCSU Legends Car Model Mass and Inertia Properties

Parameter Description

Value

Chassis (Sprung Mass) Mass Chassis Inertia Tensor Left Front Spindle (Unsprung) Mass Left Front Spindle Inertia Tensor Left Front Wheel Inertia Tensor Right Front Spindle (Unsprung) Mass Right Front Spindle Inertia Tensor Right Front Wheel Inertia Tensor Rear Axle Assembly (Unsprung) Mass Rear Axle Inertia Tensor Left Rear Wheel Inertia Tensor Right Rear Wheel Inertia Tensor

23.3 slugs {{ 110, 0, 0 },{ 0, 775, 0 },{ 0, 0, 775 }} slug-ft2 2.07 slugs {{ 2, 0, 0 },{ 0, 2, 0 },{ 0, 0, 2 }} slug-ft2 {{ 4.5, 0, 0 },{ 0, 8.97, 0 },{ 0, 0, 4.5 }} slug-ft2 2.07 slugs {{ 2, 0, 0 },{ 0, 2, 0 },{ 0, 0, 2 }} slug-ft2 {{ 4.5, 0, 0 },{ 0, 8.97, 0 },{ 0, 0, 4.5 }} slug-ft2 6.74 slugs {{ 30, 0, 0 },{ 0, 10, 0 },{ 0, 0, 30 }} slug-ft2 {{ 4.5, 0, 0 },{ 0, 8.97, 0 },{ 0, 0, 4.5 }} slug-ft2 {{ 4.5, 0, 0 },{ 0, 8.97, 0 },{ 0, 0, 4.5 }} slug-ft2

the vehicle with respect to the center of gravity of each mass, several reference points on the chassis, rear suspension and spindles were chosen to ease the measurement process. The measurements were made with respect to these reference points and then converted to the body fixed coordinate systems. The chassis reference points on the NCSU Legends car were taken to be the bottom edge of the inside corners of the frame rails which make up the boundaries of the passenger compartment floor. These are roughly the same points used to check chassis height during the vehicle setup process. The measurements were taken by placing the vehicle on jack stands above a large sheet of paper. A plum bob was used to find the projection of key points onto the ground plane. A combination square was used to determine the height of the joints above the ground plane. The measurements for the front spindles were taken with the spindles still attached to the vehicle but with the springs and the shocks removed from the vehicle. The wheels were supported so that their position relative to the chassis was roughly the same as when the vehicle was sitting on the ground. This was done primarily to reduce

154

Table 9.5 - Suspension Spring and Damper Properties

Parameter Description Front Suspension Spring Free Length Spring Preload Length Spring Rate (Linea r) Spring Rate (Cubic) Damping Coeff. (Jounce, Linear) Damping Coeff. (Jounce, Cubic) Damping Coeff. (Rebound, Linear) Damping Coeff. (Rebound, Cubic) Rear Suspension Spring Free Length Spring Preload Length Spring Rate (Linear) Spring Rate (Cubic) Damping Coeff. (Jounce, Line ar) Damping Coeff. (Jounce, Cubic) Damping Coeff. (Rebound, Linear) Damping Coeff. (Rebound, Cubic)

Left Side

Right Side

0.833 ft 0.24 ft 2340 lb/ft 0 lb/ft 3 150 lb-s/ft 0 lb(s/ft)3 150 lb-s/ft 0 lb(s/ft)3

0.833 ft 0.35 ft 2460 lb/ft 0 lb/ft 3 150 lb-s/ft 0 lb(s/ft)3 150 lb-s/ft 0 lb(s/ft)3

0.833 ft 0.462 ft 2400 lb/ft 0 lb/ft 3 150 lb-s/ft 0 lb(s/ft)3 150 lb-s/ft 0 lb(s/ft)3

0.833 ft 0.433 ft 2100 lb/ft 0 lb/ft 3 150 lb-s/ft 0 lb(s/ft)3 150 lb-s/ft 0 lb(s/ft)3

measurement error by aligning the spindle so that it was roughly parallel to the ground, thereby making the ground plane parallel to the y-axis of the spindle coordinate system. The rear axle was supported in the same manner. The geometric vehicle data is presented in Table 9.1, Table 9.2 and Table 9.3. The mass properties of the various bodies were determined as follows. The rear axle was weighed by removing the springs and shocks and weighing each wheel with the chassis supported on jack stands. The front wheel and spindle assemblies were weighed in a similar manner. Inertias for the wheels were estimated using geometric data from the wheels and tires and using the appropriate material densities (see Appendix B for a more detailed explanation). The inertia of the sprung mass was estimated by extrapolating the data found in (Garrot, 1988) to the appropriate total vehicle weight. The mass and inertia

155

properties for the various bodies which comprise the vehicle model are presented in Table 9.4. Spring rates were assumed to be linear and were computed based on the physical dimensions of the springs (Shigley, 1989). Damping rates were estimated using measured data from Winston Cup car shock absorbers which were scaled by the ratio of the vehicle masses. This scaling can be justified if one models the vehicle as a simple one degree of freedom mass-spring-damper system. The damping ratio for such a system is written as

ζ =

C 2 Mω n

(9.2.1)

Equating the damping ratios for the Winston Cup car and the Legends Car and assuming that the natural frequencies are the same gives

ζ =

CWC CLC = 2 MWCω n 2 M LCω n

(9.2.2)

or M  CLC =  LC CWC  MWC 

(9.2.3)

Tire Data The tire data consists of the geometric parameters which describe the tire as well as data describing the tractive properties of the tire. The tire data used for simulating the NCSU Legends car came from several sources. The geometric data was obtained by direct measurement of the tire. The tire equivalent spring stiffness was estimated based on values which appear in the literature. Damping in the tire was neglected. The relaxation length values and low speed damping threshold were set based on the recommendations of 156

Table 9.6 - Miscellaneous Tire Model Parameters: Geometric Data, Slip Equation Parameters and Normal Force Parameters.

Parameter Name

Parameter Description

Value

Fz0 R0 K0 K1 C0 C1 RLX RLY DCX DCY DAMPVEL

nominal wheel load tire radius (no load) radial tire stiffness (linear coefficient) radial tire stiffness (cubic coefficient) radial tire damping (linear coefficient) radial tire damping (cubic coefficient) longitudinal relaxation length lateral relaxation length longitudinal slip low speed damping coefficient lateral slip angle low speed damping coefficient low speed damping threshold

350 lbs 0.938 ft 18000 lb/ft 0 lb/ft 3 0 lb-s/ft 0 lb(s/ft)3 0.3 ft 3.0 ft 0.8 0.8 0.5 ft/s

Bernard and Clover (Bernard, 1995). Table 9.6 summarizes the values used for these parameters. The lateral and longitudinal force generation characteristics were modeled using tire data taken from an information packet supplied by BFGoodrich for use by the collegiate Legends car racing teams. The data packet is reproduced in Appendix C. The tires are designated as BFGoodrich Comp TA HR4 “Legends Edition” and are derived from a passenger car tire design. Tire data was provided for three tire pressures (15 psi, 25 psi and 35 psi) and four wheel loads (1125 lbs, 900 lbs, 675 lbs and 450 lbs). The tire data packets consists of the following plots: Cornering Stiffness versus Load, Lateral Force versus Slip Angle and Aligning Moment versus Slip Angle. No data was provided for longitudinal force versus longitudinal slip and no data was provided for combined slip operation. Under typical racing conditions the wheel loads for the NCSU Legends car are in the 150 lb to 750 lb range. Thus, the most useful data sets are those that were taken at 675 lbs and 450 lbs normal load. The 450 lb normal load data curves were used since they

157

Table 9.7 - Delft ’97 Tire Model Parameters: Pure Longitudinal Slip Equation

Parameter Name

Parameter Description

Value

P_CX1 P_DX1 P_DX2 P_EX1 P_EX2 P_EX3 P_EX4 P_HX1 P_HX2 P_KX1 P_KX2 P_KX3 P_VX1 P_VX2

longitudinal force curve shape factor longitudinal coefficient of friction longitudinal coefficient of friction - normal load dependence longitudinal force curvature factor longitudinal force curvature factor - normal load dependence longitudinal force curvature factor - normal load dependence longitudinal force curvature factor - longitudinal slip asymmetry longitudinal force horizontal offset longitudinal force horizontal offset - normal load dependence initial longitudinal force stiffness initial longitudinal force stiffness - normal load dependence initial longitudinal force stiffness - normal load dependence longitudinal force curve vertical offset longitudinal force curve vertical offset - normal load dependence

1.5 1.2 0 0 0 0 0 0 0 20 1 0 0 0

are most representative of the load conditions seen by the tire. No attempt was made to model the variation of the tire model parameters with normal load. The Magic Formula Tire Model, discussed in Chapter 7, consists of several types of curves which are fitted to the available tire data. One of the advantages of using the Magic Formula Tire Model is that the constants which appear in the equations have physical significance and can be easily obtained from tire performance plots, such as the ones in Appendix C. Due to the lack of more detailed data on the Legends car tires a number of the model parameters were set to zero. These parameters are primarily responsible for modeling the more subtle dependencies of the tire forces on such things as normal load variation and tire camber angles. The parameters which are responsible for modeling the slight asymmetries between the positive and negative slip regions were also set to zero.

158

Table 9.8 - Delft ’97 Tire Model Parameters: Pure Lateral Slip Equation

Parameter Name

Parameter Description

Value

P_CY1 P_DY1 P_DY2 P_DY3 P_EY1 P_EY2 P_EY3 P_EY4 P_HY1 P_HY2 P_HY3 P_KY1 P_KY2 P_KY3 P_VY1 P_VY2 P_VY3 P_VY4

lateral force curve shape factor lateral coefficient of friction lateral coefficient of friction - normal load dependence lateral coefficient of friction - camber angle dependence lateral force curvature factor lateral force curvature factor - normal load dependence lateral force curvature factor - camber angle dependence lateral force curvature factor - camber angle dependence lateral force horizontal offset lateral force horizontal offset - normal load dependence lateral force horizontal offset - camber angle dependence initial cornering stiffness initial cornering stiffness - normal load dependence initial cornering stiffness - camber angle dependence lateral force curve vertical offset lateral force curve vertical offset - normal load dependence lateral force curve vertical offset - camber angle dependence lateral force curve vertical offset - camber and normal load dep.

1.207 0.95 0 0 -0.932 0 0 0 0 0 0 20.05 1 0 0 0 0 0

The parameters for the longitudinal pure slip curve and the lateral pure slip curve are shown in Table 9.7 and Table 9.8 respectively. Since the steering system model isn’t sensitive to aligning torque, no effort was made to model the aligning torque and all of the associated coefficients are set to zero. The overturning moment and the rolling resistance were also ignored and their associated coefficients were set to zero. The coefficients in the preceding tables describe the lateral and longitudinal force generation characteristics of the tires under pure slip conditions (i.e. acceleration/braking or cornering, but not both at the same time). It is also necessary to characterize the relationship between lateral and longitudinal force generation under conditions of combined slip (simultaneous acceleration/braking and cornering). Due to the lack of data on combined slip tire properties, the general shapes for the combined slip curves were estimated based on plots obtained from the literature for other tires [e.g. Gillespie, 1992].

159

The parameters for the combined slip equations are shown in Table 9.9.

Track Data The road surface data provides the geometric description of the road surface. The data for the road model used for the simulation results presented in this chapter was obtained from measurements made on the Kenley, NC race track. The measurements were made by pacing off the track dimensions. The bank angles were measured using a protractor and a bubble level. Although the dimensions of the runoff apron were measured, and are included in the drawing below, the runoff apron was not included in the final road surface model since the car doesn’t normally drive it. A widened version of the track model, utilizing a 100 ft track width instead of the measured 65 ft width, was created and used during optimization to make it easier for the optimizer to find valid steering profiles.

Table 9.9 - Delft ’97 Tire Model Parameters: Combined Slip Equations

Parameter Name

Parameter Description

Value

R_BX1 R_BX2 R_CX1 R_HX1 R_BY1 R_BY2 R_BY3 R_CY1 R_HY1 R_VY1 R_VY2 R_VY3 R_VY4 R_VY5 R_VY6

longitudinal force - longitudinal slip dependence longitudinal force - longitudinal slip dependence longitudinal force - minimum value coefficient longitudinal force - horizontal offset lateral force - lateral slip dependence lateral force - lateral slip dependence lateral force - lateral slip dependence lateral force - minimum value coefficient lateral force - horizontal offset lateral force - vertical offset lateral force - vertical offset, normal load dependence lateral force - vertical offset, camber dependence lateral force - vertical offset, lateral slip dependence lateral force - vertical offset, longitudinal slip dependence lateral force - vertical offset, longitudinal slip dependence

1.0 0.5 9 0 16.5 0 0 1.04 0 0 0 0 0 0 0

160

Figure 9.1 - The Schematic of the Kenley, NC Race Track

Driver Model Data The driver data consists of the gains and other parameters which describe the driver’s response to command inputs. The driver model parameters for the simulation were set as shown in Figure 9.1. The minimum and maximum steering angle parameters specify the range of optimization for the steering profile points. The optimization range is asymmetric about zero because all of the turns on the road course are left turns. The minimum and maximum velocity parameters specify the optimization range for the velocity profile points. The delay constant parameter (DM_VEL_T_DELAY) is the time constant for a first order ODE which is used to filter the steering input function before it is applied to the vehicle. The filter acts to remove the high frequency components of the steering input function which keeps the integration algorithm from reducing it’s step size

161

unnecessarily. The roles of the remaining constants were discussed in the chapter on driver modeling.

9.3 Model Chassis Setup Once the vehicle measurements have been made it is usually necessary to make fine adjustments to the suspension in the same way that a real race car is set up prior to a race. The setup process includes setting the toe angle, camber angle and caster angle for the front wheels, setting the track width for the front suspension, squaring the rear axle with respect to the chassis, setting the left side and right side wheel base, setting the cross weight and finally, setting the frame heights. The vehicle setup process for a real car is performed in several steps. It is frequently necessary to repeat the sequence of steps because an adjustment made to one part of the car can upset an adjustment made elsewhere. The cross weight adjustments are very sensitive to irregularities in the surface, thus, all adjustments should be made with the car sitting on a flat and level surface. The first step is to set the ride height to the desired value. This is done by adjusting the position of the spring seat on each of the four coil-over shocks used on the NCSU Legends car. This effectively lengthens or shortens the shock body. Lengthening the shock body raises the associated corner of the car. To set the ride height, all four shocks are adjusted until the desired frame heights are achieved at each of the four corners of the car. The next step is to align the rear axle. This is done by adjusting the trailing links to square the rear axle with respect to the chassis and by changing the length of the panhard

162

Table 9.10 - Driver Model Parameters

Parameter Name

Parameter Description

Value

DM_STEER_MAXANGLE DM_STEER_MINANGLE DM_VEL_MAXSPEED DM_VEL_MINSPEED DM_VEL_T_DELAY DM_VEL_T_PREVIEW DM_VEL_LAMBDA DM_VEL_M_EFF DM_VEL_BRAKEBIAS DM_VEL_MAXACCEL DM_VEL_MAXDECEL DM_VEL_STTC_OFFSET DM_VEL_STTC_GAIN

maximum allowed steering angle for steering profile opt. minimum allowed steering angle for steering profile opt. maximum allowed speed for velocity profile optimization minimum allowed speed for velocity profile optimization time constant for high frequency filter for steering angle preview (look ahead) time velocity error feedback loop gain factor command acceleration gain (effective vehicle mass) front/rear brake bias constant maximum allowed acceleration maximum allowed deceleration traction control threshold traction control gain factor

1.0472 rad -0.5236 rad 140 ft/s 40 ft/s 0.05 sec 0.003 sec 2.0 70 slugs 0.37 100 ft/s2 -30 ft/s2 0.8 0.25

rod to set the axle’s lateral position with respect to the chassis. It is also possible, and sometimes desirable, to apply a small amount of steering angle to the rear axle by adjusting the relative lengths of the trailing links. The rear wheel track width can be increased by using shims between the wheel mounting flange and the wheel rim but is otherwise fixed. The next step is to align and to position the front wheels relative to the rear axle and the chassis. The camber angle, the caster angle and the toe angle for each of the front wheels is set by adjusting the lengths of the control arms and strut rods. The wheel base can be adjusted by lengthening or shortening the three trailing links supporting the rear axle or by adjusting the lengths of the strut rods supporting the front spindles. The track width is set by increasing or decreasing the lengths of the control arms.

163

Once the basic vehicle setup geometry has been achieved, the cross weights can be set and the frame heights can be rechecked. To set the cross weight, the four wheels of the vehicle are placed on scales and the lengths of the shocks on opposite corners of the car are adjusted together to achieve the desired weight distribution. It may be necessary to iterate between adjusting the frame height and the cross weight. If the suspension geometry changes significantly due to raising or lowering the car it may also be necessary to reset the suspension geometry. The setup process for the model follows a nearly identical sequence of steps. The simulation code used to generate the results discussed in this chapter also has the capability of solving for the equilibrium position of the system. By placing the model on a flat road surface and finding the equilibrium position, the wheel weights and the front and rear suspension alignment parameters can be determined. By iteratively adjusting the suspension link lengths, in the same manner as discussed in the preceding paragraphs, the desired vehicle setup may be achieved. The vehicle setup parameters used for the NCSU Legends car simulation are shown in Table 9.11. The suspension link lengths presented in

Table 9.11 - Vehicle Setup Parameters

Parameter Description

Value

Cross Weight Left Front and Right Front Frame Heights Left Rear and Right Rear Frame Heights Left Front Camber Right Front Camber Front toe angle Front Track Width Rear Track Width Left Side Wheel base / Right Side Wheel base Rear Axle Steer Angle

50.5% 3.625” 3.875” 3.0 degrees (out at top) -5.0 degrees (in at top) 0.20” out 50.9” 49.9” 73” ± ¼” 0.0 degrees

164

the suspension geometry tables presented earlier in this chapter are the result of performing this setup process on the vehicle model.

9.4 Simulation Results Prior to starting a vehicle optimization it is necessary to establish a baseline for comparison purposes. The baseline run is identical to the vehicle optimization run in that it optimizes the steering profile and the velocity profile to obtain the best possible lap time and the minimum path tracking error. The difference is that the baseline optimization run does not include vehicle design parameters as degrees of freedom. The optimizer setup, the steering profile configuration and results, the velocity profile configuration and results, and the simulation results are discussed in the following sections.

Optimizer and Cost Function Computation Setup The optimizer parameters were set as shown in Table 9.12. The cost function value is computed as CF = WPosition CFPosition + WVelocity CFVelocity + WLap _ Time CFLap _ Time

(9.4.4)

where the CFi represent the cost function components and the Wi are weighting

Table 9.12 - Optimization and Cost Function Computation Parameters

Parameter Description

Value

Iffco_RMaxH Iffco_StartH Iffco_RMinH Iffco_Fscale Iffco_Restarts CFW_PathErr CFW_VelErr CFW_LapTime CFW_LapPenalty

0.5 0.5 2.0e-006 100 4 1.0 0.001 1.0 40.0

165

parameters. The cost function for the baseline optimization run utilized all three cost function components (position, velocity and lap time) using the weights shown in the table. The cost function scale parameter was set to a value of 100 which is approximately equal to the weighted sum of the maximum values of the cost function components: CFMAX = 10 . (66.8)+ 0.001(100.0) + 10 . (40.0) ≈100

(9.4.5)

The maximum values for the cost function components were obtained as follows. The maximum lateral position error is equal to the maximum possible lateral position penalty as discussed in section 9.1. The maximum velocity error is a conservative estimate based on experience with the simulation code. The maximum lap time cost is specified by the user via the CFW_LapPenalty parameter (see Table 9.12) and represents a time to be assigned to laps which are not completed. The velocity error cost function component, which is equal to the time integral of the instantaneous error between the vehicle velocity and the desired velocity, is weighted to minimize it’s contribution to the overall cost function. Experience has shown that using a larger weight for this cost function component can prevent the optimizer from improving the lap time. This is because increases in the peak speeds can lead to increased velocity errors when the car is nearing it’s acceleration limits and/or it’s tractive limits. Setting the weight to zero would be fine except that the simulation code uses the cost function component weight as a flag to determine which parameters (in this case the speed profile parameters) to include in the optimization. Thus, it is necessary to use a small, but nonzero, value for the cost function weight.

166

For a reasonable well optimized steering profile, for the Kenley, NC track, the lap time cost function component is roughly one order of magnitude larger than the path error cost function component. The weights for both of these components are set to unity. This has the tendency to emphasize lap time in the cost function minimization process. Experience has shown that preferentially weighting the lap time is necessary in order to achieve the fastest possible laps. It is believed that this is a result of the increase in the path error cost function component, which normally occurs when changes are made to the velocity profile, which has the tendency to discourage the optimizer from changing the velocity profile.

Optimal Steering Profile Configuration The optimized steering profile for the Kenley, NC track, which is shown in Figure 9.2, consists of 20 steering angle points. Fourteen of these points are fixed with respect to the track so that their position is not included in the optimization process. The remaining six points can be repositioned by the optimization algorithm as necessary. There are a total of 26 optimization degrees of freedom (2 DOF for each of the 6 unconstrained points plus 1 DOF for the remaining 14 constrained points) for this steering profile. The large number of optimization degrees of freedom makes the solution process very slow. To alleviate this problem, the first optimization run was made using a version of the steering profile which took advantage of the periodicity of the steering profile. The second half of the steering profile was made identical to the first half, with the exception of being offset to the second half of the track. This reduced the number of optimization parameters to 13 and enabled a solution to be found relatively rapidly. This solution was 167

45.0 40.0 35.0

Angle (degrees)

30.0 25.0 20.0 15.0 10.0 5.0 0.0 -5.0 -10.0 0

200

400

600

800

1000

1200

1400

1600

1800

2000

Track Position (ft) Steering Wheel Angle

Unconstrained SW Profile Point

Constrained SW Profile Point

Figure 9.2 - The Steering Profile for the Kenley, NC Simulation

used as an initial starting point for the full 26 parameter aperiodic steering profile optimization run.

Optimal Velocity Profile Configuration The optimized velocity profile for the Kenley, NC track, which is shown in Figure 9.3, consists of four points. The symmetry of the track allows for periodicity of the velocity profile. Accordingly, the second pair of velocity profile points is based on a shifted image of the first pair of points. This reduces the number of optimization parameters from a total of eight parameters (speed and position for each of the four points) to four parameters. Unlike the steering profile, the position for all of the velocity profile points is free to be optimized. This allows the optimizer to pick the best point at

168

which to begin braking for the corner and the best point at which to begin accelerating out of the corner. The velocity profile used as a starting point for this optimization was the result of numerous prior simulation and optimization runs and, as such, is somewhat atypical of a normal initial velocity profile. In the absence of prior knowledge of the vehicle’s performance limits, one starts with a velocity profile with moderate speeds which are guaranteed to allow the car to complete a lap. Including the lap time term in the overall cost function allows the optimizer to gradually increase the speeds at which the car navigates the road course until it is no longer possible for the car to remain on the desired path.

75.0

70.0

Velocity (mph)

65.0

60.0

55.0

50.0

45.0 0

200

400

600

800

1000

1200

1400

1600

Track Position (ft) Desired Long Vel

Aperiodic Points

Periodic Points

Figure 9.3 - The Velocity Profile for the Kenley, NC Simulation

169

1800

2000

Speed Control Algorithm Performance The driver speed control performance is quite good, as shown in Figure 9.4, which compares the prescribed velocity profile and the actual vehicle velocity profile. The best lap time obtained for the baseline simulation was 21.4 seconds. This is slower than the lap times obtained for the real vehicle on the Kenley, NC track. The best times obtained to date by NCSU Legends car team members are in the 17.4 second to 18.0 second range. The discrepancy in the lap times is not due to the inability of the model to match the speeds attained by the real vehicle. The maximum velocity on the back stretch and the minimum velocity in the corner were measured for the real Legends car using a radar gun. The minimum speed was approximately 40 mph to 44 mph, depending on the driver. The

75

70

Velocity (mph)

65

60

55

50

45 0.0

5.0

10.0

15.0

20.0

Time (sec) Vehicle Velocity

Prescribed Velocity

Figure 9.4 - A Comparison of the Prescribed Velocity and the Actual Vehicle Velocity

170

20.0

Acceleration (ft/s^2)

15.0

10.0

5.0

0.0

-5.0 0.0

5.0

10.0 Time (sec)

15.0

20.0

Figure 9.5 - The Vertical Acceleration of the Sprung Mass (Sprung Mass Coordinate System)

maximum speed was between 71 mph and 78 mph. The optimized speed profile for the simulation, shown in Figure 9.4, matches the measured speed ranges quite well. It is possible that the sawtooth pattern used for the prescribed velocity profile does not represent the velocity profile of the real vehicle very well. The real vehicle may be spending a greater amount of time at the high speeds, thus reducing the lap time. Another possible explanation for the discrepancy is that the combined slip characteristics of the tire model may not be correct. The combined slip behavior was modeled in an ad hoc fashion due to the lack of tire data for the combined slip loading. The real Legends car may exhibit greater cornering power, allowing it to exit the corner at greater speeds, again reducing lap times.

171

0.30

0.25

Slip (unitless)

0.20

0.15

0.10

0.05

0.00

-0.05

20.0

15.0

10.0

5.0

0.0

-0.10

Time (sec) LF Long Slip

RF Long. Slip

LR Long. Slip

RR Long. Slip

Figure 9.6 - The Longitudinal Wheel Slip Percentages

One final possible explanation is that the track dimensions could be in error. As mentioned earlier, the track was measured by pacing off it’s dimensions. An over-estimate of 10%, which would not be out of the question given the methods used, could easily cause the lap times to be 2 seconds slower. The measured bank angles at the midpoints of the corners have been confirmed to be within half a degree. The bank angles at the corner entrance, corner exit, and at the midpoint of the straight away have not been confirmed and could be error. Aside from the slow lap times, the performance of the driver speed controller is very good. There are two deviations from the prescribed velocity profile which occur at

172

Figure 9.7 - The Vehicle Position at 9.0 seconds (Exiting Turn 2)

9.0 seconds and at 19.5 seconds. These deviations are a result of the traction control algorithm reducing the drive torque in order to limit wheel spin. The longitudinal wheel slip percentages are shown in Figure 9.6. The traction control threshold was set to activate at a longitudinal slip of 0.2. The rear wheel slip percentages exceed the traction control threshold value at 9.0 seconds and at 19.5 seconds which matches the deviations from the velocity profile. The reason for the loss of traction can be seen in Figure 9.5 which contains a plot of the vertical acceleration of the vehicle. The vertical acceleration takes on negative values at the times in question. This is due to a slight crest in the track surface which occurs at the exit of each of the turns. A snapshot of the vehicle’s position at 9.0 seconds

173

800

Normal Load (lb)

700 600 500 400 300 200 100 0.0

5.0

10.0

15.0

20.0

Time (sec) Left Front

Right Front

Left Rear

Right Rear

Figure 9.8 - The Tire Normal Loads

is shown in Figure 9.7. Note that a widened version of the Kenley, NC track is shown in the figure. The use of a widened track makes it easier for the optimizer to find an optimal steering profile. The location of the driver path is, in actuality, very close to the inside of the curve when viewed on the unmodified track. The plot of the normal loads on each tire (Figure 9.8) shows similar evidence of a lack of normal force leading to the loss of traction.

Steering Control Performance The performance of the steering controller can be analyzed by looking at the lateral position error of the vehicle with respect to the prescribed driver path over the course of the simulation. This plot is shown in

Figure 9.9. The overall performance of the driver

controller is quite good with an average tracking error of approximately 2.5 ft. The worst case tracking error is approximately 9 ft and occurs at approximately 18 sec. The reason

174

for this large position error becomes apparent upon inspection of the related vehicle performance plots. The most interesting of these plots are the yaw velocity plot, the lateral acceleration plot, and the steering wheel angle plot. The yaw velocity plot shows peaks at approximately 6.0 seconds and at approximately 16.5 seconds. The car is in the middle section of the turns from 4 seconds to 8 seconds (turns 1 and 2) and from 15 seconds to 19 seconds (turns 3 and 4). The yaw velocity peaks are significantly higher than the average yaw velocity through the turn which indicates that there is a loss of traction. When viewing an animation of the vehicle’s motion, it is possible to see the rear end of the car slip toward the outside of the corner at the times noted above.

10.0 9.0

Lateral Position Error (ft)

8.0 7.0 6.0 5.0 4.0 3.0 2.0 1.0 0.0 0.0 Average Position Error = 2.84 ft

5.0

10.0

15.0

Time (sec)

Figure 9.9 - The Vehicle Lateral Position Error

175

20.0

40 35

Angular Velocity (deg/s)

30 25 20 15 10 5 0 -5 0.0

5.0

10.0

15.0

20.0

Time (sec)

Figure 9.10 - The Yaw Velocity

Inspection of the lateral acceleration plot ( Figure 9.11) provides addition confirmation of the loss of traction at the rear wheels. The lateral acceleration plot shows two strong peaks at 6.0 seconds and at 17.0 seconds. These peaks have magnitudes of approximately 1.25 G’s which exceeds the tractive capabilities of the tires by a significant margin. The steering angle plot ( Figure 9.11) again confirms the loss of traction. The plot shows two minima, located in the center of the turns, at 6.2 seconds and at 16.9 seconds. The steering wheel angle drops to approximately 4 degrees at both minima. The steering angle at the beginning and end of each of the turns is approximately 42 degrees. The minima in the steering profile occur because the driver is counter-steering the car to

176

45 40

Steering Wheel Angle (deg) or Lateral Acceleration (ft/s^2)

35 30 25 20 15 10 5 0 -5 -10 0.0

5.0

10.0

15.0

20.0

Time (sec) Lateral Acceleration

Steering Wheel Angle

Unconstrained SW Profile Points

Locked SW Profile Points

Figure 9.11 - The Steering Wheel Angle and the Lateral Acceleration

compensate for the loss of traction at the rear wheels in order to remain on the desired driver path.

9.5 Vehicle Optimization Results The baseline optimization results generated above were used as a starting point for a relatively simple vehicle setup optimization. The left rear spring rate and the left rear spring free length were selected as the optimization parameters. Modification of these parameters affects the vehicle’s cross weight and the vehicle’s weight transfer

Table 9.13 - Vehicle Suspension Parameter Optimization Ranges

Parameter Name

Parameter Description

Lower Bound

Upper Bound

DP_RS_LSP_FL DP_RS_LSP_K0

left rear spring free length left rear spring rate

0 inches 166.7 lbs/in

14.4 inches 250.0 lbs/in

177

characteristics, and thus, affects the oversteer/understeer behavior of the car during cornering. The optimization was run using the parameter set which was used for the baseline optimization. The only change was the addition of the two vehicle suspension parameters mentioned above. The optimization ranges for the suspension parameters are shown in Table 9.13. The initial values and optimized values for the suspension parameters are shown in Table 9.14. As can be seen from the optimization results, the changes were minimal. Although the initial setup is known to be a good setup for the real vehicle when running on the real track, it is unlikely that this setup is the best possible setup. Also, given the lack of accurate tire data and other unavoidable discrepancies between the model and the real vehicle, there are most likely some handling differences between the model and the real vehicle which would lead to differences in the optimal configuration. Also of note is the fact that the car demonstrated a tendency to oversteer in the apex of the corners during the baseline optimization as noted in the preceding section. Based on these observations, it seems clear that the vehicle optimization failed to improve the vehicle setup. This result is not entirely unexpected given the difficulty of simultaneously optimizing the steering profile and the velocity profile which has been encountered in the

Table 9.14 - Vehicle Suspension Parameter Optimization Results

Parameter Name

Parameter Description

Original Value

Optimized Value

DP_RS_LSP_FL DP_RS_LSP_K0

left rear spring free length left rear spring rate

9.996 inches 200 lbs/in

9.997 inches 200.06 lbs/in

178

past and which was noted in the preceding section. The explanation for this behavior is that making significant changes to the vehicle setup (or to the prescribed velocity profile) causes the vehicle’s path to deviate from the desired path by a significant amount. If the steering profile has been reasonably well optimized already, this leads to an increase in the path error cost function which discourages the optimizer from making changes to the vehicle setup unless it is able to simultaneously correct the steering profile to eliminate the driver path tracking error. This problem, and possible solutions, are discussed in the following section.

9.6 Recommendations for Future Research The vehicle model itself seems to be capable of reproducing the fundamental behavior of the NCSU Legends car on the Kenley, NC track. At this point no detailed experimental validation of the model has been performed. It is recommended that the car be fitted with a data acquisition system so that comparisons can be made to the model predictions. The data acquisition system should include sensors for recording the same type of data from the real car as was presented in the preceding sections for the model car. This list of sensors includes the following devices: steering wheel angle position sensor, three axis accelerometer, three axis angular velocity sensor, wheel velocity sensors (all four wheels), and finally an engine RPM sensor. Another useful sensor would be a high sample rate high accuracy global positioning system or another system with equivalent functionality. This would allow determination of the driver path around the track as well as providing a better estimate of

179

vehicle velocity which, in combination with the wheel velocities, could be used to compute the longitudinal slip ratios for the wheels. By driving the vehicle around the inner periphery and outer periphery of the track at low speeds it could also be used for measuring the track itself. At present, no engine model is included in the simulation. The current arrangement provides whatever power is necessary to achieve the desired acceleration, up to the limits of adhesion of the tires. While it is not an unreasonable simplification when applied to the Legends car, it would not be terribly difficult to limit acceleration based on the available engine power at the current engine RPM (this can be calculated based on rear wheel speed, gear ratio, etc.). Another area of possible impro vement deals with the braking system. The current algorithm simply limits the maximum deceleration of the vehicle to avoid lockup. Applying a more intelligent control algorithm, perhaps one similar to the currently implemented traction control algorithm, could improve braking performance entering the corners which would allow the optimizer to maintain straight away speed a bit longer or to increase the maximum speed. The braking system on the real car includes a device which causes the brake pressure at the rear wheels to lag the brake pressure applied to the front wheels. This can assist in preventing rear wheel lockup when the brakes are applied rapidly. This aspect of the braking system is not modeled in the simulation. It would also be useful to include brake bias as a parameter in the optimization to extract the maximum possible performance from the braking system.

180

The biggest improvements can be made in the driver control algorithms. The use of the optimizer to handle the steering and throttle control tasks, in addition to optimizing the vehicle design parameters, interferes with the optimization process and can prevent the code from improving the performance of the vehicle. The large number of parameters being optimized with this approach severely degrades the performance of the optimizer. It is highly desirable to separate the driver control task from the vehicle parameter optimization. One possible approach would be to break the road course into a fairly fine mesh of segments which can be traversed sequentially using a shooting-method approach: Given that the vehicle is on course at the beginning of the current segment, determine the steering input required to arrive at the beginning of the next segment. An interpolating polynomial of fairly low order could be used to generate the steering profile over the length of the segment. The steering angle at the beginning of the segment is known (from traversing the preceding segment). The steering input at the end of the segment, and perhaps one or two other parameters, could be treated as the unknowns in the shooting problem. The velocity at which the vehicle traverses the segment could be handled in a similar manner. This approach has the advantage of guaranteeing that the car will either compete the lap in an acceptable manner or that it will fail to complete the lap. Either way, the task of driving the vehicle has been removed from the main optimization loop. The only potential problem with this approach is that the vehicle response typically lags the steering input by a significant amount (see Figure 9.11 - The Steering Wheel Angle and the Lateral Acceleration for an example of this). When the length of the segments gets too small the

181

steering inputs from the current segment may have a strong effect on the handling of the vehicle when it is traversing the succeeding segment. This could lead to more and more over-compensation as the vehicle crosses from one segment to the next and result in instability of the driver controller. It may be possible to improve the stability to some extent by enforcing a path tangency criterion in addition to a path offset criterion as the vehicle traverses the segment.

182

Bibliography Allen, R. Wade and McRuer, Duane T. (1979). “The Man/Machine Control Interface Pursuit Control”, Automatica, 15:683-686. Allen, R. Wade. (1982). “Stability and performance analysis of automobile driver steering control”. SAE Technical Paper Series , Paper #820303. Allen, R. Wade, Henry T. Szostak, Theodore J. Rosenthal and Donald E. Johnston. (1986). “Test Methods and Computer Modeling for the Analysis of Ground Vehicle Handling”, SAE Technical Paper Series , Paper #861115. Allen, R. Wade, Theodore J. Rosenthal, Henry T. Szostak. (1987a). “Steady State and Transient Analysis of Ground Vehicle Handling”, SAE Technical Paper Series , Paper #870495. Allen, R. Wade, Henry T. Szostak, and Theodore J. Rosenthal. (1987b). “Analysis and computer simulation of driver/vehicle interaction”, SAE Technical Paper Series , Paper #871086. Allen, R. Wade, Theodore J. Rosenthal, Henry T. Szostak. (1988). “Analytical Modeling of Driver Response in Crash Avoidance Maneuvering Volume I: Technical Background”, U.S. Department of Transportation, NHTSA, DOT HS 807 270. Allen, R. Wade, Henry T. Szostak, Theodore J. Rosenthal, and David H. Klyde. (1990). “Field Testing and Computer Simulation Analysis of Ground Vehicle Dynamic Stability”, SAE Technical Paper Series , Paper #900127. Allen, R. Wade, Henry T. Szostak, Theodore J. Rosenthal, David H. Klyde, and K. J. Owens. (1991). “Characteristics Influencing Ground Vehicle Lateral/Directional Stability”, SAE Technical Paper Series , Paper #910234. Allen, R. Wade, Theodore J. Rosenthal, David H. Klyde, K. J. Owens and Henry T. Szostak. (1992). “Validation of Ground Vehicle Computer Simulations Developed for Dynamics Stability Analysis”, SAE Technical Paper Series , Paper #920054. Allen, R. Wade, Theodore J. Rosenthal. (1993). “A Computer Simulation Analysis of Safety Critical Maneuvers for Assessing Ground Vehicle Dynamic Stability”, SAE Technical Paper Series , Paper #930760.

183

Allen, R. Wade, Theodore J. Rosenthal, Jeffrey R. Hogue. (1996). “Modeling and Simulation of Driver/Vehicle Interaction”, SAE Technical Paper Series , Paper #960177. Antoun, R. J., P. B. Hackert, M. C. O’Leary, and A. Sitchin. (1986). “Vehicle Dynamic Handling Computer Simulation - Model Development, Correlation, and Application Using ADAMS”, SAE Technical Paper Series , Paper #860574. Bakker, E., L. Nyborg and H.B. Pacejka. (1987). “Tyre Modeling for Use in Vehicle Dynamics Studies”, SAE Paper #870421. Bakker, E., H.B. Pacejka and L. Lidner. (1989). “A New Tire Model with an Application in Vehicle Dynamics Studies”, SAE Paper #890087. Becker, G. , H. Fromm and H. Maruhn. (1931) “Schwingungen in Automobillenkungen” (“Vibrations of the Steering Systems of Automobiles”), Krayn, Berlin. Bekey, G. A., G. O. Burnham and J. Seo. (1977) “Control Theoretic Models of Human Drivers in Car Following”, Human Factors, 19(4):399-413. Bell, Steven C., W. Riley Garrott, John R. Ellis, Y. C. Liao. (1987). “Suspension Testing Using the Suspension Parameter Measurement Device”, SAE Technical Paper Series , Paper #870577. Bernard, James E. (1973). “Some time-saving methods for the digital simulation of highway vehicles”. Simulation, 21(6):161-165, December 1973. Bernade, James E. and C. L. Clover. (1994). “Validation of Computer Simulations of Vehicle Dynamics”, SAE Technical Paper Series , Paper #940231. Bernard, James E. and C.L. Clover. (1995). “Tire Modeling for Low-Speed and HighSpeed Calculations”, SAE Paper #950311. Bryan, G. H. (1911). “Stability in Aviation”, Macmillan & Co., London. Broulheit, G. (1925). “La Suspension de la Direction de la Voiture Automobile: Shimmy et Danadinement’ (“The Suspension of the Automobile Steering Mechanism: Shimmy and Tramp”), Société des Ingénieurs Civils de France, bulletin 78, 1925. Bundorf, R. Thomas. (1967). “The influence of vehicle design parameters on characteristic speed and understeer”. Transactions of the Society of Automotive Engineers, SAE 670078. Clover, Chris L. and James E. Bernard. (1993). “The Influence of Lateral Load Transfer Distribution on Directional Response”, SAE Technical Paper Series , Paper #930763.

184

Crossman, E. R. F. W., H. Szostak and T. L Cesa. (1966). “Steering Performance of Automobile Drivers in Real and Contact-Analog Simulated Tasks”, Human Factors Society Tenth Annual Meeting , Oct. 1966. Chrstos, Jeffrey P. (1991). “A Simplified Method for the Measurement of Composite Suspension Parameters”, SAE Technical Paper Series , Paper #910232. Dickson, J. G. and A. J. Yardley. (1993). “Development and Application of a Functional Model to Vehicle Development”, SAE Technical Paper Series , Paper #930835. Donges, Edmund. (1978). “A Two-Level Model of Driver Steering Behavior”, Human Factors, Vol. 20(6), pp. 691-707. Etheridge, Mark C. (1998). “Preliminary Performance of Carbon-Carbon Valves in High Speed Pushrod Type Valve Trains”, Master’s Thesis, North Carolina State University, Raleigh, NC. Evans, R. D. (1935). “Properties of tires affecting riding, steering, and handling”. Journal of the Society of Automotive Engineers, 36(2):41. Garrot, W. Riley, Douglas L. Wilson and Richard A. Scott. (1981). “Digital Simulation For Automobile Maneuvers”, Simulation , pp. 83-91, September 1981. Garrot, W. Riley, Michael W. Monk and Jeffrey P. Chrstos. (1988). “Vehicle Inertial Parameters - Measured Values and Approximations”, SAE Technical Paper Series , Paper #881767. Genta, Giancarlo. (1997). “Motor Vehicle Dynamics: Modeling and Simulation”, ISBN 9810229119. Gillespie, Thomas D. (1992). “Fundamentals of Vehicle Dynamics”, Society of Automotive Engineers, ISBN 1-56091-199-9. Goland, Martin and Frederick Jindra. (1961). “Car Handling Characteristics”. Automobile Engineer , 51(8):296-302, August 1961. Gordon, D. A. (1966a). “Experimental Isolation of Drivers’ Visual Input”, Public Roads, 33(12):266-273. Gordon, D. A. (1966b). “Perceptual Basis of Vehicular Guidance”, Public Roads, 34(3):53-68. Gruening, James and James E. Bernard. (1996). “Verification of Vehicle Parameters for Use in Computer Simulation”, SAE Technical Paper Series , Paper #960176. Guo, K. and H. Guan. (1993). “Modelling of Driver/Vehicle Directional Control System”, Vehicle System Dynamics , 22:141-184.

185

Heydinger, Gary J., W. Riley Garrot, Jeffrey P. Chrstos, and Dennis A. Guenther. (1990). “A Methodology for Validating Vehicle Dynamics Simulations”, SAE Technical Paper Series, Paper #900128. Heydinger, Gary J., Paul A. Grygier, and Seewoo Lee. (1993). “Pulse Testing Techniques Applied to Vehicle Handling Dynamics”, SAE Technical Paper Series , Paper #930828. Heydinger, Gary J., Nicholas J. Durisek, David A. Coovert Sr., Dennis A. Guenther and S. Jay Novak. (1995). “The Design of a Vehicle Inertia Measurement Facility”, SAE Technical Paper Series , Paper #950309. Hickner, G. B., J.G. Elliot and G.A. Cornell. (1971). “Hybrid Computer Simulation of the Dynamic Response of a Vehicle With Four-Wheel Adaptive Brakes”, SAE Technical Paper Series, Paper #710225. Hoffman, E. R., (1975). “Human Control of Road Vehicles”, Vehicle System Dynamics , 5(1-2):105-126. Huang, Feng, J. Roger Chen and Lung-Wen Tsai. (1993). “The Use of Random Steer Test Data for Vehicle Parameter Estimation”, SAE Technical Paper Series , Paper #930830. Jindra, Frederick. (1976). “Mathematical Model of Four-Wheeled Vehicle for Hybrid Computer Vehicle Handling Program”, Deparment of Transportation - National Highway Traffic Safety Administration, DOT HS-801-800, January 1976. Kim, Dojoong. (1990). “Dynamics and Optimal Design of High Speed Automotive Valve Trains”, Doctoral Thesis, North Carolina State University, Raleigh, NC. Kondo, M. and A. Ajimine. (1968). “Drivers Sight Point and Dynamics of the DriverVehicle System Related to it.”, SAE Technical Paper Series , Paper #680104. Kortüm W. and W. Schiehlen. (1985). “General Purpose Vehicle System Dynamics Software Based on Multibody Formalisms”, Vehicle System Dynamics , 14:229-263. Kortüm W. and R. S. Sharp. (1993). “Multibody Computer Codes in Vehicle System Dynamics”, Supplement to Vehicle System Dynamics , Volume 22. Lin, Y. and W. Kortüm. (1991). “Identification of System Physical Parameters for Vehicle Systems with Nonlinear Components”, Vehicle System Dynamics , 20:354-365. Lu, Zhengyu, Andrzej G. Nalecz, Kenneth L. d’Entremont. (1993). “Development of Vehicle-Terrain Impact Model for Vehicle Dynamics Simulation”, SAE Technical Paper Series, Paper #930833. McHenry, R. and N. Deleys. (1968). “Vehicle Dynamics in Single Vehicle Accidents”, Technical Report CAL No. VJ-2251-V-3, Cornell Aeronautical Laboratory Inc., December 1968.

186

McRuer, Duane T., R. Wade Allen, David H. Weir and Richard H. Klein. (1977). “New Results in Driver Steering Control Models”, Human Factors, 19(4):381-397. Milliken, William F. and David W. Whitcomb. (1956). “General introduction to a programme of dynamic research”. Proceedings of the Automobile Division; The Institution of Mechanical Engineers , (7):287-309. Mimuro, Tetsushi, Takahiro Maemura and Hiroshi Fujii. (1993). “Development and Application of the Road Profile Measuring System”, SAE Technical Paper Series , Paper #930257. Modjtahedzadeh, A. and R. A. Hess. (1993). “A model of driver steering control behavior for use in assessing vehicle handling qualities”. Transactions of the ASME: Journal of Dynamic Systems, Measurement, and Control, 115:456-464, September 1993. Mori, Yoshinori, Hironobu Matsushita, et al. (1991). “A Simulation System for Vehicle Dynamics Control”, SAE Technical Paper Series , Paper #910240. Morman, Kenneth N. (1977). “Non-Linear Model Formulation for the Static and Dynamic Analyses of Front Suspensions”, SAE Technical Paper Series , Paper #770052. Mousseau, C. W., M. W. Sayers and D. J. Fagan. (1991). “Symbolic Quasi-Static and Dynamic Analyses of Complex Automobile Models.”, Proceedings, 12 th International Association for Vehicle System Dynamics Symposium on the Dynamics of Vehicles on Roads and Tracks, Lyon, France, Aug. 1991. Murphy, Ray W. (1970). “A Hybrid Computer System for the Simulation of Vehicle Dynamics”, SAE Technical Paper Series , Paper #700154. Nalecz, Andrzej G. (1987). “Investigation into the Effects of Suspension Design on Stability of Light Vehicles”, SAE Technical Paper Series , Paper #870497. Nalecz, Andrzej G. et al. (1988). “Effects of Light Truck and Roadside Characteristics on Rollover”, Final Report for First Period, NHTSA - US DOT Contract No. DTNH22-89C-07005, August 23, 1988. Nalecz, Andrzej G. (1992). “Development and Validation of Light Vehicle Dynamics Simulation (LVDS)”, SAE Technical Paper Series , Paper #920056. Nikravesh, Parviz E. (1988). “Computer-Aided Analysis of Mechanical Systems”, ISBN 0-13-164220-0 025, Prentice-Hall, Inc., Englewood Cliffs, New Jersey. Oberdieck W., B. Richter and P. Zimmerman. (1979). “Adaptation of Mathematical Vehicles Models to Experimental Results”, Proceedings of the 6 th IAVSD Symposium on the Dynamics of Vehicles on Roads and Tracks, Berlin, Germany, pp. 367-378, September 3-7, 1979.

187

Okada, T., T. Takiguchi, M. Nishioka, and G. Utsumomiya. (1973). “Evaluation of Vehicle Handling and Stability by Computer Simulation at the First Stage of Vehicle Planning”, SAE Technical Paper Series , Paper #730525. Olley, M. (1934). “Stable and Unstable Steering”, Unpublished Report, General Motors, 1934. Olley, M. (1937). “Suspension and Handling”, Engineering Report, General Motors Corp., May 1937. Pacejka, H. B. and I. J. M Besselink. (1997). “Magic Formulate Tyre Model with Transient Properties”, Vehicle System Dynamics Supplement , Vol. 27, pp. 234-249. Pacejka, H. B. and E. Bakker. (1992). “The Magic Formula Tyre Model”, Tyre Models for Vehicle Dynamics Analysis , 1st International Colloquium on Tyre Models for Vehicle Dynamics Analysis, Supplement to Vehicle System Dynamics , pp. 1-18. Vol. 21. Petersen, Michael R. and John M. Starkey. (1996). “Nonlinear Vehicle Performance Simulation with Test Correlation and Sensitivity Analysis”, SAE Technical Paper Series , Paper #960521. Radt, H. S. and W. F. Milliken Jr. (1960a). “Motions of Skidding Automobiles”, SAE Summer Meeting, Paper #205A, June 1960. Radt, H. S. and W. F. Milliken Jr. (1960b). “Exactly What Happens When an Automobile Skids?”, SAE Journal, pp. 27-33, December 1960. Radt, H. S. (1964). “Analysis of Nonlinear Problems in Vehicle Handling”, Advanced Dynamic Institute, Wayne State University, September 1964. Roland, R. D. and T. B. Sheridan. (1966). “A Normative Model for Control of Vehicle Trajectory in an Emergency Manoeuvre”, Paper presented at Annual Meeting of the Highway Research Board, Jan. 1966. Roland, R. D. and T. B. Sheridan. (1967). “Simulation Study of the Driver’s Control of Sudden Changes in Previewed Path”, Mass. Inst. Tech., Department of Mechanical Engineering, Report DSR 74920-1, 1967. Sayers, Michael W. and Paul S. Fancher. (1993). “Hierarchy of Symbolic ComputerGenerated Real-Time Vehicle Dynamics Models”, Transportation Research Record, 1403:88-97. Segel, Leonard. (1956a). “Research in the Fundamentals of Automobile Control and Stability”, SAE National Summer Meeting, Atlantic City, June 5, 1956.

188

Segel, Leonard. (1956b). “Theoretical prediction and experimental substantiation of the response of the automobile to steering control”. Proceedings of the Automobile Division; The Institution of Mechanical Engineers , (7):310-330, 1956. Segel, Leonard. (1966). “On the Lateral Stability and Control of the Automobile as Influences by the Dynamics of the Steering System”, Transactions of the ASME: Journal of Engineering for Industry , August 1966. Sharp, R. S. (1994). “The Application of Multi-body Computer Codes to Road Vehicle Dynamics Modeling Problems”, Proceedings of the Institution of Mechanical Engineers, Part D, Journal of Automobile Engineering , 208(D1):55-61. Shigley, Joseph Edward and Charles R. Mischke. (1989). “Mechanical Engineering Design - 5th Edition”, McGraw-Hill, Inc. Smiley, R. F. and W. B. Horne. (1958). “Mechanical Properties of Pneumatic Tires with Special Reference to Modern Aircraft Tires”, NACA, TN 4110, 1958 or NASA, TR R-64, 1960. Speckhart, Frank H. (1973). “A Computer Simulation for Three-Dimensional Vehicle Dynamics”, SAE Technical Paper Series , Paper #730526. Society of Automotive Engineers. (1965). “Vehicle Dynamics Terminology SAE J670a”. An updated version is available (SAE J670e). Tiffany, N. O., G. A. Cornell and R. L. Code. (1970). “A Hybrid Simulation of Vehicle Dynamics and Subsystems”, SAE Technical Paper Series , Paper #700155. Walker, G. E. L. (1950). “Directional stability, a study of factors involved in private car design”. Automobile Engineer , 40(530,533):281-370. Weir, David H., C. P. Shortwell, and W. A. Johnson. (1968a). “Dynamics of the automobile related to driver control”. Transactions of the Society of Automotive Engineers , SAE 680194. Weir, David H. and Duane T. McRuer. (1968b). “A Theory for Driver Steering Control of Motor Vehicles”, Highway Research Record, 247(2):7-27. Weir, David H. and Richard J. DiMarco. (1978). “Correlation and Evaluation of Driver/Vehicle Directional Handling Data”, SAE Technical Paper Series , Paper #780010. Whitcomb, David W. and William F. Milliken. (1956). “Design implications of a general theory of automobile stability and control”. Proceeding of the Automobile Division; The Institution of Mechanical Engineers , (7):367-391. Yoshimoto, K. (1969). “Simulation of Man-Automobile Systems by the Driver’s Steering Model with Predictability”, Bulletin of the J.S.M.E. , 12(51):495-500.

189

Appendix A Useful Derivatives In this section the derivatives of commonly appearing quantities are calculated. Derivatives which are zero are not listed explicitly. Derivatives are taken with respect to each generalized coordinate, with respect to each generalized velocity and with respect to time for most quantities.

Angular Velocity Derivatives Derivatives are calculated in this section for the angular velocity

R

ω R/E of a

rotating coordinate system R with respect to the inertial coordinate system E. The angular velocity of a rotating coordinate system with respect to the inertial coordinate system is

R

ω R /E

− β1 r&  = 2 [L]β = 2 − β2  − β3

β0 − β3 β2

β3 β0 − β1

& β − β2   &0   β  β1  &1  β2 β0  &    β3 

(A.1)

& − ββ & − ββ & &  β β 2 3 + β3β2  0 &1 1 &0 & − ββ & = 2 β0β2 + β1β3 − β2β 0 3 1 β β & + ββ & &  0 &3 − β1β 2 2 1 − β3β0  The only nonzero derivatives are those taken with respect to the βi and the β&i .

∂ R ω R /E ∂β0

& β & β 0 1 0 0  &0  1     β1  & = 2 0 0 1 0 &  = 2 β2  β2 β    &3  0 0 0 1  &   β3 

190

(A.2)

∂ R ω R /E ∂β1

− 1  = 2 0  0

∂ R ω R /E ∂β2

0  = 2 − 1  0

∂ R ω R /E ∂β3

0 0  = 2 0 − 1  − 1 0

∂ ω R /E ∂β&0

− β1  = 2 − β2  − β3

∂ ω R /E ∂β&1

− β1  = 2 − β2  − β3

∂ R ω R /E ∂β&2

− β1  = 2 − β2  − β3

∂R ω R /E ∂β&3

− β1  = 2 − β2  − β3

R

R

0 0 0

0 0 1

0 0 −1

&  β 0 &  − β    0 0 &    β1  & 1  &  = 2  β3  β2 − β  0  &2   &   β3 

(A.3)

0 0 0

&  β & − β   − 1 &0  3    β1  & 0  &  = 2 − β0  β2 β  0  &1   &   β3 

(A.4)

1 0 0

& β 0 &  β    0 &  &2   β1  0  &  = 2 − β1  β2 − β  0  &0   &   β3 

(A.5)

β0

β3

− β3

β0

β2

− β1

β0

β3

− β3

β0

β2

− β1

β0

β3

− β3

β0

β2

− β1

β0

β3

− β3

β0

β2

− β1

1 − β1  − β2       0 β1    = 2 − β2  0    β0  − β3   0

(A.6)

0  β0  − β2       1 β1    = 2 − β3  0    β0   β2   0

(A.7)

0  β3  − β2       0 β1    = 2  β0  1    β0  − β1   0

(A.8)

0 − β2  − β2       0 β1    = 2  β1  0    β0   β0   1 

(A.9)

Note that the time derivative, taken with respect to the inertial coordinate system, is

191

E

∂R ω R /E d ∂R ω R /E  R d ∂R ω R /E  R     ω = + × R /E dt  ∂β&i  dt  ∂β&i  ∂β&i ∂R ω R /E  R ∂R ω R /E    ω R /E × =− + ∂βi ∂β&i  

(A.10)

where R

∂ R ω R /E d ∂ R ω R /E    & = − dt  ∂β ∂βi i

(A.11)

The time derivative of the angular velocity can be found using the vector differentiation theorem E

R

R d R d R ω R /E ) = ( ( ω R /E ) + R ω R /E ×R ω R /E dt dt R d R = ( ω R /E ) dt

(A.12)

r d R &r& &r& & ]β& = 2 [L] β ω R /E ) = 2 [L]β + 2 [L ( dt

R

ω& R /E

− β1  = 2 − β2  − β3

β0

β3

− β3

β0

β2

− β1

&&  β && − β β && && &&  (A.13) β0β   − β2 &&0  1 0 − β2β3 + β3β2  1   β1  && + β β && && && β1  &&  = 2β0β 2 1 3 − β2 β0 − β3β1  β2 β β && + β β && &&  β0   0 &&3 − β1β 2 2 1 − β3β0   &&   β  3

Transformation Matrix Derivatives Derivatives of the coordinate transformation matrix [E C R ] , which transforms a vector from its representation in the rotating system to its representation in the inertial system, are calculated in this section. The coordinate transformation matrix is expressed in terms of Euler Parameters. The definition of the transformation matrix is given below in terms of two linear transformation matrices. The resulting nonlinear transformation matrix depends only on the βιand time (indirectly).

192

β20 + β12 − 12 β1β2 − β0β3 β1β3 + β0β2    [E C R ] = [G ][L]T = 2 β1β2 + β0β3 β20 + β22 − 12 β2β3 − β0β1  2 2 1  β1β3 − β0β2 β2β3 + β0β1 β0 + β3 − 2   − β1  [G ] = − β2  − β3

β0

− β3

β3 − β2

β0 β1

− β1 β2    − β1  ; [L] = − β2  β0   − β3

β0

β3

− β3 β2

β0 − β1

− β2   β1  β0  

(A.14)

(A.15)

The derivatives of the transformation matrix with respect to the βι are 2β0 ∂[E C R ]  = 2  β3 ∂β0  − β2

− β3

 0 ∂[E C R ]  = 2  β1 ∂β2  − β0

β1

2β0 β1 2β2 β3

β2   − β1  2β0  

2β1 β2 ∂[E C R ]  = 2 β2 0 ∂β1  β3 β0

β3   − β0  0  

β0   β3  0 

0 ∂[ E C R ]  = 2 β0 ∂β3  β1

β1   β2  2β3  

− β0 0 β2

(A.16)

The first time derivative of the transformation matrix is d & ] = [G & ][L]T + [G ][L & ]T = 2[G & ][L]T = 2[G ][L & ]T [ E C R ] =[ E C R dt & − β β&  β0β&0 + β1β&1 − β2β&2 − β3β&3 − β0β&3 + β1β&2 + β2β 1 3 0  & β0β&0 − β1β&1 + β2β&2 − β3β&3 = 2 β0β3 + β1β&2 + β2β&1 + β3β&0 & & & & & & & & − β β  0 2 + β1β3 − β2 β0 + β3β1 β0 β1 + β1β0 + β2 β3 + β3β2

β0β&2 + β1β&3 + β2β&0 + β3β&1   − β0β&1 − β1β&0 + β2 β&3 + β3β&2  &  β0β&0 − β1β&1 − β2 β&2 + β3β 3  (A.17)

Note that ~ ] & ] =[ C ][R ω [E C R E R R /E

(A.18)

where [R ω~ R/ E ] is the skew symmetric matrix associated with the angular velocity vector R ω R/ E . R ω R/ E is the angular velocity of the non-inertial coordinate system with respect to the inertial coordinate system expressed in terms of the unit vectors of the noninertial coordinate system. The result of the product of the skew symmetric matrix ~ ] with a vector x is a vector containing the cross product of R ω and x. [R ω R/ E R/ E

193

& ] with respect to the β are The derivatives of [E C R i

& & β − β 0 3 & ] ∂[E C  R & & β0 = 2  β3 ∂β0 & & − β  2 β1 & β & − β 2 1 & ] ∂[E C  R & & = 2  β1 β2 ∂β2 & & − β  0 β3

&  β 2  & − β1  &  β 0  & β 0 & β3 & − β

    2

& β 1 & ] ∂[ E C  R & = 2 β2 ∂β1 & β 3 & − β 3 & ] ∂[E C  R & = 2  β0 ∂β3 & β  1

& β 2 & − β1 & β 0 − −

& β 0 & β3 & β 2

&  β 3  & − β0  & − β 1 & β 1  & β2  &  β 3

(A.19)

& are & ] with respect to the β The derivatives of [ E C R i

 β0 & ] ∂[E C  R = 2  β3 & ∂β0  − β2

− β3 β0 β1

− β2 β1 & ] ∂[E C  R = 2  β1 β2 & ∂β2  − β0 β3

β2   − β1  β0   β0   β3  − β2  

β1 & ] ∂[E C  R = 2 β2 & ∂β1  β3 − β3 & ] ∂[E C  R = 2  β0 & ∂β3   β1

β2 − β1 β0 − β0 − β3 β2

β3   − β0  − β1   β1   β2  β3  

(A.20)

Note that & ] ∂[ C & ] d ∂[ E C R E R = dt ∂β&i ∂βi

194

(A.21)

The second time derivative of the transformation matrix is d2 & ][L & ]T + 2[G &&][L]T = 2[G & ][L &]T + 2[G ][L &&]T [ C R ] = 2[G 2 E d t &2 + β &2 − β &2 − β&2 − 2β& β& + 2β& β& β 2β&0β&2 + 2β&1β&3  0 1 2 3 0 3 1 2 &β &  &β & &2 − β&2 + β&2 − β&2 − 2β& β& + 2β && ] = 2 2β& β& + 2β β [E C R 1 2 0 1 2 3 0 1 2 3   0 3 &β & && &β & & & &2 − β &2 − β&2 + β&2  − 2β 2β β 0 2 + 2β1β3 0 1 + 2β2β3 0 1 2 3  && + β β && && && && + β β && && &&  β0β − β0β 0 1 1 − β2β2 − β3β3 3 1 2 + β2β1 − β3β0  && && && && && − β β && && && β0β + 2 β0β 3 + β1β2 + β2β1 + β3β0 0 1 1 + β2β2 − β3β3 && + β β && && && β β && && && && − β0β 2 1 3 − β2β0 + β3β1 0 1 + β1β0 + β2β3 + β3β2 

(A.22)

&& + β β && && &&  β0β 2 1 3 + β2β0 + β3β1 && − β β && && &&  − β0β 1 1 0 + β2β3 + β3β2  && − β β && && &&  β0β 0 1 1 − β2β2 + β3β3  Alternatively, the second time derivative of the transformation matrix may be found by differentiating equation A.18

(

)

d2 d & d ~ ] [ C R ] = [E C [ E C R ][R ω R]= R/ E 2 E d t dt dt ~ ] + [ C ][ R ω ~& ] & ][R ω =[ E C R R/ E E R R/E ~ ~ ~& =[ C ][R ω ][R ω ] + [ C ][R ω E

R

R/ E

R/ E

E

R

(A.23) R/ E

]

and substituting the time derivative of the angular velocity [R ω~& R/ E ]

[ R

0  ~  & & & & & && − β β && ω R /E = 2 β0β3 − β1β2 + β2β 1 3 0  && − β β && && && − β0β 2 1 3 + β2β0 + β3β1 

]

&& + β β && && && − β0β 3 1 2 − β2β1 + β3β0 0 & & & & && + β β && β0β1 − β1β0 − β2β 3 3 2

&& + β β && && &&  β0β 2 1 3 − β2β0 − β3β1 && + β β && && &&  − β0β 1 1 0 + β2β3 − β3β2   0 

to obtain the same result as before.

195

(A.24)

Appendix B Wheel Inertia Estimate The inertia of the wheel brake and tire assemblies can be estimated by approximating the assembly as a collection of cylindrical disks and cylindrical shells of varying thicknesses and materials. The inertias for the basic shapes are computed below.

Thin Cylindrical Disk The general formula for computing the inertia of a body about a particular axis a is I a = ∫ra dm 2

(B.1)

V

where ra is the distance of the mass element dm from the axis of rotation. The inertias for objects possessing cylindrical symmetry are most easily computed in cylindrical coordinates. The mass element in this case is given by

dm = ρ r dr dθ dz

(B.2)

where ρ is the density of the material. For a thin disk, with thickness T, outside radius Ro and inside radius Ri, there are only three unique inertias. They are as follows:

I xy =

T / 2 2 π Ro

∫ ∫∫r

2

cos(θ)sin (θ)ρ r drdθ dz = 0

(B.3)

− T / 2 0 Ri

I xx = I yy =

T / 2 2π Ro

∫ ∫∫r

2

cos2 (θ)ρ r drdθ dz =

− T / 2 0 Ri

I zz =

T / 2 2π Ro

∫ ∫∫r

− T / 2 0 Ri

2

ρ r drdθ dz =

π π ρT Ro 4 − Ri 4 = ρT Do 4 − Di 4 4 64

(

)

(

π π ρT Ro 4 − Ri 4 = ρT Do 4 − Di 4 2 32

(

196

)

(

)

)

(B.4)

(B.5)

The z-axis is the axis of rotational symmetry. The x and y axes are in the plane of the disk and are orthogonal to each other and to the z-axis.

Thin Walled Cylindrical Shell The inertias for a thin-walled cylindrical shell are computed in the same manner as for the disk. The shell is assumed to have a mean radius R, thickness T and a length L. The z-axis is assumed to be coincident with the axis of rotational symmetry. The x and y axes are constructed to be orthogonal to one another and to the z-axis. The volume of the cylindrical shell is 2 2  T  T  V = ∫ ∫ ∫r drdθ dz = π LR +  − R −   = 2π RLT 2  2   − L / 2 0 R− T / 2   L / 2 2π R + T / 2

(B.6)

and thus, the total mass is

m = 2πρ RLT

(B.7)

The inertias are

I xx = I yy =

L / 2 2π R + T / 2

∫ ∫ ∫ρ( z

2

)

+ r 2 sin 2 (θ) r drdθ dz

− L / 2 0 R− T / 2

π ≈ ρ L 3 RT + πρLR 3T 6 m m = 6 R 2 + L2 = 3D 2 + 2 L2 12 24

(

)

(

(B.8)

)

4 4 π  T  T  I zz = ∫ ∫ ∫r ρ r drdθ dz = LρR +  − R −   2  2  2  − L / 2 0 R− T / 2   π = Lρ 4 R3T + RT 3 ≈2πρ LTR3 2 m = mR 2 = D 2 4 L / 2 2π R + T / 2 2

[

]

197

(B.9)

Rotating Assembly Model The rotating assembly can be modeled as a collection of thin cylindrical shells and thin disks. The rotating assembly includes the brake disk or brake drum, the wheel mounting flange assembly, the wheel and the tire. The wheel center is modeled as a thin disk and the wheel rim is modeled as a thin cylindrical shell. The portion of the brake disk swept by the brake pads is modeled as a thin disk. The brake disk mounting flange is also modeled as a thin disk and the material connecting the brake disk to the mounting flange is modeled as a thin walled cylinder. The tire sidewalls are modeled as thin disks and the tread surface is modeled as a thin walled cylinder. The wheel mounting flange is modeled as a thin disk. Note that, when computing the inertia of the complete assembly, it is necessary to apply the parallel axis theorem to obtain the proper result. The densities for the materials which make up the rotating assembly are easily obtained with the exception of the tire density. The tire density can be estimated by measuring it’s weight and computing an approximate volume. This approach ignores the uneven distribution of mass in the tire due to the presence of the steel belts but can provide a reasonable estimate. If a tire is not available for measurement it is possible to estimate the weight of the tire by measuring the weight of the wheel and tire together and subtracting the weight of the wheel rim. The weight of the wheel rim is computed by estimating it’s volume and multiplying by the appropriate material density.

198

Appendix C Tire Data

199

200

201

202

203

204