Tire and force distribution modeling and validation for wheel loader applications
John Spencer and Bernhard Wullt
Mechanical Engineering, master's level 2017
Luleå University of Technology
Department of Engineering Sciences and Mathematics
Acknowledgements
We would like to thank our supervisors, Andre Fernandez and David Berggren, at Volvo CE as well as Jan-Olov Aidanp¨a¨a at LTU. Your help and support during the thesis work has been of high importance for our work. We would also like to thank Lennart Skogh and his test department, who helped us out a lot with the measurements and the testing. Last but not least, we would like to give a big thank you to Kausihan Selvam and Auayporn Elfving at Volvo CE for giving so freely of their time. You have been a great help to us! Thanks!
Abstract
This thesis describes the development of a machine force distribution estimator and the calibration of tire models for wheel loader applications.
Forces generated in the contact patch between the tires and the ground are crucial for under- standing and controlling machine dynamics. When it’s not possible to directly measure these contact patch forces they are estimated from other sensor data. Validated models of the contact patch forces are also used in machine dynamics simulations and are very relevant to model based development. Vehicle dynamics is of crucial importance to the automotive industry. In contrast the modeling of these forces has not been very important to the construction equipment industry and as such wheel loaders haven’t been studied as much as conventional cars.
In order to model the forces in the contact patch, suitable tire models have been studied and calibrated. The tire models have been calibrated by using two different sources, field data and test rig data. Two steady state tire models were chosen for the field data. These were the brush model and the Magic Formula. The resulting fit from the data was not ideal, but this was due to that the given data were of low accuracy. The Magic Formula was used for the test rig data, which gave a good overall fit. The results from the test rig measurements were then used in a transient model, the single contact point model, and the MF-Tyre software. The models were implemented in Simulink and were validated against experimental data. They showed good correspondence, but deviated for some levels of slip.
Another important aspect of the wheel loader is the force distribution over the entire machine.
Two estimators have been developed, one to estimate the vertical forces on each tire, the normal force estimator and one to calculate the turning behavior due to different force outputs on the tires, the turning torque estimator. The normal force on each tire is information that is important for the tire model, but it can also be used to estimate when the wheel loader risk tipping on its side. The turning torque estimation is useful for control systems to optimize the driving behavior of the machine.
Compared against measured data from an actual wheel loader the normal force estimator showed a high accuracy in estimating the individual wheel vertical forces. The turning torque estimator could estimate the behavior of the torque but had problems when estimating the magnitude.
Table of Content
1 Introduction 1
1.1 Thesis work . . . 1
1.1.1 Scope . . . 2
1.2 Literature review . . . 3
2 Theory 4 2.1 Tire model: Tire fundamentals . . . 4
2.1.1 Introduction to tire quantaties . . . 4
2.1.2 Tire forces and moments . . . 5
2.2 Tire model: Tire models . . . 7
2.2.1 The Magic Formula . . . 7
2.2.2 MFSWIFT . . . 9
2.2.3 The brush model . . . 11
2.2.4 Single contact point model . . . 13
2.3 Tire model: Algorithms for parameter estimation . . . 13
2.3.1 Differential evolution . . . 14
2.4 Force distribution: Normal force . . . 14
2.4.1 Defining the basis for the normal force distribution model . . . 14
2.4.2 Finding the center of gravity . . . 15
2.4.3 Shifts in center of gravity . . . 16
2.4.4 Normal force on each tire . . . 20
2.4.5 Roll-over . . . 25
2.4.6 Pile Entry . . . 25
2.5 Force distribution: Acceleration from wheel sensors . . . 26
2.6 Force distribution: Turning torque . . . 27
3 Method 30 3.1 Tire model: Test vehicle . . . 30
3.2 Field measurments . . . 30
3.2.1 Tire used for measurements . . . 30
3.2.2 Measurement systems . . . 30
3.2.3 Experiments . . . 32
3.2.4 Complementary field test data . . . 32
3.3 Tire model: Test rig measurements . . . 33
3.3.1 Tire used for measurements . . . 33
3.3.2 Main setup . . . 33
3.3.3 Pure longitudinal force characteristics . . . 35
3.3.4 Relaxation length . . . 36
3.3.5 Validation . . . 36
3.4 Tire model: Optimization routine for parameter identification of tire models . . . 36
3.5 Tire model: Parameter and model identification of field measurements . . . 36
3.5.1 Reference tire model . . . 37
3.5.2 Magic Formula parameters . . . 37
3.5.3 Brush model . . . 38
3.6 Tire model: Data processing of test rig measurements . . . 38
3.6.1 Processing of force characteristics measurements . . . 38
3.6.2 Processing of relaxation length measurements . . . 38
3.7 Tire model: Parameter and model identification of test rig measurements . . . 39
3.8 Tire model: Model validation of test rig measurements . . . 39
4 Results 41 4.1 Tire model: Processed field measurements . . . 41
4.2 Tire model: Parameter and model identification from field measurements . . . 41
4.2.1 Magic Formula parameters . . . 41
4.2.2 Brush model . . . 44
4.2.3 Comparison of tire models . . . 46
4.3 Tire model: Parameter and model identification from test rig measurements . . . 47
4.3.1 Magic Formula parameters for nominal load and inflation pressure . . . 47
4.3.2 Magic Formula parameters for variation of load and nominal inflation pres- sure . . . 48
4.3.3 Comparison of resulting fitted curves for variation of loads . . . 50
4.3.4 Magic Formula parameters for nominal load and variation of inflation pres- sure . . . 51
4.3.5 Comparison of resulting fitted curves for variation of inflation pressures . . 52
4.3.6 Relaxation length . . . 53
4.4 Tire model: Model validation of models from test rig measurements . . . 53
4.5 Force distribution: Normal force estimation results . . . 54
4.5.1 Normal force estimator vs Simulation . . . 54
4.5.2 Transient torque case . . . 54
4.5.3 Transient steering angle case . . . 55
4.5.4 Using the old measurement data data . . . 56
4.5.5 Using newer complimentary data . . . 58
4.5.6 Roll over . . . 60
4.5.7 Pile entry . . . 61
4.6 Force distribution: Turning torque results . . . 62
5 Discussion 64 5.1 Tire model: Field measurements . . . 64
5.2 Tire model: Test rig measurements . . . 64
5.3 Force distribution: Normal force . . . 65
5.3.1 Compared against simulation . . . 65
5.3.2 Compared against the old tests . . . 65
5.3.3 Compared against the new tests . . . 65
5.3.4 Roll over . . . 66
5.3.5 Pile entry . . . 66
5.4 Force distribution: Turning torque . . . 66
6 Conclusions 67 6.1 Tire model: Field measurements . . . 67
6.2 Tire model: Test rig measurements . . . 67
6.3 Force distribution: Normal force estimation . . . 68
6.4 Force distribution: Turning torque . . . 68
7 Future work 69 7.1 Tire model . . . 69
7.2 Force distribution estimator . . . 69
References 72
Appendices 73
A Tire model: Data processing for field measurements 73
B Tire model: Processed field measurement data 79
C Tire model: Structural properties from field measurements 83
D Tire model: Parameters and bounds for field measurements 84
E Tire model: Parameters and bounds for test rig measurements 85
F Kalman filter 86
Variable list
Name Symbol Unit
Half contact patch length ac m
Longitudinal distance between wheel and COG a m Lateral distance between wheel and COG b m
Stiffness factor B -
Longitudinal stiffness cpx N/m2
Shape factor C -
Longitudinal slip stiffness CF κ N
Lateral slip stiffness CF α N
Normalized change in load dfz -
Normalized change in inflation pressure dpi -
Peak factor D N
Slope factor E -
Unit vector e -
Longitudinal Force Fx N
Steady state force Fxss N
Lateral Force Fy N
Vertical Fore Fz N
Nominal load Fzo N
Normal force from asphalt data Fza N
Normal force from gravel data Fzg N
Gravity constant g m/s2
Weighting function G -
Population for differential evolution Gi -
Cost function J -
Stiffness Kxκ N
Moment around the longitudinal axle Mx Nm
Moment around the lateral axle My Nm
Moment around the vertical axle Mz Nm
Mass m kg
Curbed mass, L-60 mc kg
Bucket mass, L-60 mb kg
Gross mass, L-60 mg kg
Added mass, L-60 ml kg
Mass of tire mt kg
Revolutions per second of cardan shaft ncs rps
Inflation pressure of back tire pb bar
Pressure in piston pc bar
Inflation pressure of front tire pf bar
Inflation pressure tire pi bar
Nominal inflation pressure pio bar
Longitudinal parameters for shape factor PCx - Longitudinal parameters for peak factor PDx - Longitudinal parameters for slope factor PEx - Longitudinal parameters for stiffness factor PKx - Longitudinal parameters for inflation pressure PP x - Longitudinal parameters for horizontal offset PHx - Longitudinal parameters for vertical offset PV x -
Loaded radius r m
Unloaded radius rf m
Rotation matrix R -
Horizontal factor SH -
Vertical factor SV -
Slip point S -
Longitudinal deflection u m
Velocity v m/s
Speed of wheel V m/s
m/s Slip speed Vs m/s
Linear speed of rolling Vr m/s
Max speed for L-60 Vmax m/s
Velocity cylinder Vcyl km/h
Velocity vector V m/s
Tire contact patch width w m
Wheel position wp m
Longitudinal position x m
Longitudinal position of center of mass Xcog m
Lateral position y m
Lateral position of center of mass Ycog m Measured value in cost function Ymeasured N
Model output in cost function Ymodel N
Vertical position z m
Vertical position of center of mass Zcog m
Side-slip angle α rad
Camber angle of tire γ rad
Steering angle δ deg
Longitudinal composite function θx -
Pitch θ deg
Longitudinal slip κ -
Location of peak longitudinal force κm -
Transient longitudinal slip κ0 -
Friction coefficient, asphalt µa -
Friction coefficient, gravel µg -
Friction coefficient, concrete µc -
Peak friction µp -
Friction variable in Magic Formula µx -
Friction coefficient µ -
Radial deflection of the tire ρ m
Longitudinal theoretical slip σx -
Longitudinal theoretical slip when sliding occurs σslx -
Relaxation length σ m
Roll ψ deg
Rotational speed tire Ω rad/s
1 Introduction
Wheel loaders are used in diverse terrain types, they can be used on everything from concrete floors to muddy off-road work sites or deep inside mines running over gravel. All this while performing a range of operations from transporting huge stone blocks to stacking logs of timber. At present time Volvo Construction Equipment (VCE) has developed a functional prototype LX1, illustrated in Figure 1, which uses electric motors with individual wheel control. Conventional wheel loaders are typically slow in their dynamics due to that the drive line usually has high inertia. With individual control of each wheel the overall performance can be increased since the electric systems are de- coupled i.e. no connection between each motor. Thus each wheel can be controlled independently to achieve the maximum possible tractive force while not being affected by the ground conditions at the other wheels.
A problem for vehicles in general is unnecessary wheel spin which causes the tires to wear out faster and generates less tractive force. The tires of the wheel loader are a very expensive part, and is one of the parts that are associated with a high running cost. Thus it’s a cost that needs to be kept at a minimum. A way to control this phenomenon, is to use traction control, which has the aim to deliver the highest possible tractive force that the current condition allows. A necessary part to develop such a control algorithm is to have a simulation model that reflects the wheel loader accurately. Such a model can also be used to simulate the causes of undesirable and catastrophic effects, such as uncontrolled sliding or roll over.
A important aspect in traction control is the forces generated between the machine and the ground.
These forces can be measured directly but the available sensors for this purpose are usually not robust enough for long time use and are also too expensive. Another way to get this information, is by estimating the forces from other sensor data by using mathematical models, which is what has been done in this thesis.
The following thesis report is a combined report of two parallel works, with the goal to cali- brate a tire model and a force distribution estimator. An import part in product development is by using what’s termed as Model Based Development (MBD). This is a working process to develop functions and products through simulation, and thereby decreasing the total developing time and cost by removing the need for physical tests and prototypes. This of course puts a high demand on the accuracy of the simulation model.
Figure 1: Illustration of the LX1 prototype [1].
1.1 Thesis work
In 2015, measurements were performed where the dynamics of the wheel loader was measured in different situations by using a motion capturing system [2] together with force and moment transducers mounted on the wheels. The measurements were conducted with the primary goal to calibrate a tire model. The conducted measurement were raw field measurements and no pre- vious work at present time was known to have been performed in a similar way. Thus the used approached was an unconventional and unestablished way to perform the calibration, since its usually performed in special test rigs. However these test rigs and methods are mainly established for on-road vehicles and not for heavy off-road vehicles. Due to practical reasons, it was decided to conduct the tests as pure field tests. These measurements were used as a starting point, but further tests were conducted during the thesis work.
The normal load (vertical force) on the tires are an important input to the tire-model as it gives information regarding the available tractive forces. As the current machines are not equipped with load sensors the distribution must be estimated. This estimation must be able to handle shifts in the weight due to various dynamic phenomena. A reliable estimator of the normal load can also be used as an input for the tire model to give higher precision estimation of the forces acting on the surface. Since the wheel loader can articulate around its mid point, the normal car models used in most research papers can not be used. The center of gravity will significantly shift during turning as roughly half the machine will no longer align with the other half. Dynamic effects during a turning phase will create unbalanced weight shifts. To find the behavior of these shifts is the basis for this part of the thesis. To validate any developed model, data from wheel transducers in a com- bination of on-board systems can be used to show the real world distribution during operations.
The distribution of the the normal load on the tires will also be used to determine the likelihood of the machine tipping over. Using the normal force distribution model as a base it is also possible to calculate the force trying to turn the machine due to the variation in force on the individual tires.
At the start of the thesis work, Volvo had a simulation model of a complete wheel loader. This model was used during the thesis work as a reference tool.
1.1.1 Scope Tire model
Following list gives an overview of the different stages and goals for the tire modeling part of the thesis.
• Perform a literature review of tire kinematics and dynamics. Study tire models that has the potential to be used.
• Process the raw data from the old and newly performed experiments.
• From the processed data, perform a model and parameter identification, and develop a tran- sient and steady state tire model.
• Implement the tire models in the pre-existing simulation environment (Simulink). Perform simulations and validate the tire models.
• Implement the tire model in the real time system (ECU) of the wheel loader and evaluate from test runs.
The preconditions of the tire model are as follows.
• The model most be sufficiently fast in order to be used in the real time system.
• The tire model must be allowed to be incorporated in the pre-existing simulation environment.
• The tire model is suppose to take slip and normal force as input. The output should be the tire forces and moments.
Force distribution estimator
The work process of the force distribution estimator is:
• Literature study of weight distribution and turning moment prediction for articulated vehi- cles.
• Develop a model for finding the center of gravity for any steering angle.
• Include the dynamic shift of the gravity due to road topology and vehicle acceleration, using available sensors.
• Evaluate any special drive cases.
• Utilize the developed model to find the turning-moment from all wheels combined
• Validate the model against simulated and measured data
• Transfer the model to Simulink and run it on the plant model.
• Run the model on the on-board ECU.
• Set up plan for roll over prevention and feed-forward prediction using the steering wheel turning velocity.
Following preconditions on the force distribution estimator are:
• The normal force estimator can not use the longitudinal or lateral forces as inputs as it will be used to calculate these.
• Can not use any sensors except those for: steering angle, acceleration of vehicle, rotational wheel speed, motor torque, change in the machine articulation angle, angle of the lift arm and lift arm hydraulic cylinder pressure caused by the load
• Most be small enough in storage space to be able to be run on the wheel loaders on board computer.
1.2 Literature review
A extensive literature review was performed at the start of the thesis. In the following sections, a small but a important selection will be discussed. However, the remainder of the sources are referenced within the relevant sections within the thesis.
Tire model
The main source that was used for studying tire dynamics and models was [3]. The book gives detailed information regarding tire dynamics and different tire models. The author of the book, Hans Pacejka, is well known within the area of tire dynamics and is one the developer of the Magic Formula, which is a well known and established tire model. Due to its popularity and its abilities to describe the tire characteristics, the Magic Formula was chosen in the beginning of the thesis as the primary tire model. The book further contains useful empirical equations for the structural properties of the tire. Other books have also been used to expand the knowledge base in tire dynamics [4, 5, 6].
A big part of the thesis work has been to identify the parameters of the tire models from mea- surement data. A part of the literature review therefore dealt with finding suitable least square methods. A good article that deals with this is [7], where various algorithms were used to estimate the Magic Formula parameters. The algorithms were applied to different characteristics of the tire and the data that was used was also from different sources. The article gave a good overview and introduction to different algorithms that are suitable for the task. It also highlights some difficulties in the estimation process that were useful for the thesis work.
Force Distribution
One of the most interesting articles when performing the literature review was [8]. The authors discuss that there has been very little research in this field regarding vertical forces and roll- over. To remedy this the authors develop a model using lagrangian mechanics defined in a global coordinate system. The result is a mathematical model suitable for simulating many different drive cases. Unfortunately the proposed model requires information that will not be available for the machines that are of interest in this thesis. The general idea of how to model a wheel loader and its behavior was however used. Most of the defining equations that are used to build the model comes from [4], where the author, Reza Jazar, gives easy to understand explanations of vehicle dynamics and their effects. While the book focus on the automotive applications, the books way of building a model for how the weight of a vehicle shifts during acceleration, braking and changes in road inclination how been incorporated in the developed wheel loader model.
2 Theory
This section presents the results from the literature review conducted at the beginning of the thesis work.
2.1 Tire model: Tire fundamentals
The tire is a important component for machine dynamics. Section 2.1.1 introduces and defines some basic quantities and tire concepts. Next in section 2.1.2 the forces that act on the tire are described.
2.1.1 Introduction to tire quantaties
A tire traveling with velocity V and an angular speed Ω is presented in its general configuration in Figure 2.
Figure 2: Illustration of the tire coordinate system [4].
In [4] the tire is defined to be travelling on a flat plane termed as the ground plane. Next, the tire plane is introduced which is defined as the mid-plane corresponding to the tire [4]. A coordinate system is next defined with the origin placed at the center point of the contact patch. This co- ordinate system represents the tires local coordinate system. The x-axis is defined to be aligned with the line of intersection between the tire plane and the ground plane. The z-axis points in the direction along the normal to the ground plane. Finally, the y-axis points in the direction such that a right handed coordinate system is given. The rotation about the x, y and z-axis is termed as the roll, pitch and yaw-angle.
The following forces and moments are applied to the wheel,
• Fx, longitudinal force.
• Fy, lateral force (cornering force).
• Fz, vertical force (load).
• Mx, roll moment (overturning moment).
• My, pitch moment.
• Mz, yaw moment (aligning moment).
In [4], two different angles are presented which are the side slip and camber angle represented by α and γ respectively. The side slip angle is the angle the velocity vector, V makes with respect to the x-axis. The camber angle, γ, corresponds to the inclination of the wheel and defined as the angle the tire plane makes with the z-axis.
The loaded radius of the wheel is in [3] represented by the variable r and its unloaded radius by rf. With these radii the radial deflection, ρ, is defined [3] as
ρ = rf− r (1)
Another important radius is defined in [3], which is the effective rolling radius of the tire. It’s an important radius, and according to [3], defined as the radius which couples the longitudinal speed, Vx, and the angular speed, Ω, when the wheel is rolling freely, which means that no torque is being applied to it. From [3] it’s stated as follows
re= Vx
Ω (2)
and the radius is under free rolling conditions bounded to be between the unloaded and loaded radius of the tire [6].
r < re< rf. (3)
For the free rolling case, the radius defines the instantaneous center of rotation of the wheel and defines an slip point S which is distanced re radially from the center of the wheel [3]. This is illustrated in Figure 3.
Figure 3: Presents the slip point S and the effective rolling radius re[3].
Longitudinal and lateral slip
When a torque is applied to the wheel, the instantaneous center of rotation moves away from being located at the slip point. Relative motion is therefore given between the slip point and the wheel.
A longitudinal slip speed is given which in [3] is calculated as
Vsx= Vx− Ωre (4)
and longitudinal slip arises, which is a quantity defined by [3] as κ = −Vsx
Vx = −(Vx− reΩ)
Vx . (5)
When the wheel is driving forward, the sign of κ is positive, and positive longitudinal force, Fx, arises. If however the wheel is braked, the sign of κ turns negative and a negative longitudinal force is given. Since the previous definition of slip has the longitudinal velocity in the denominator, this gives rise to singularities when the wheel is standing still, Vx= 0, and the wheel spins, Ω 6= 0.
An alternative definition of slip [5] can instead be applied for the case of a wheel driving forward, which is as follows
κ = −Ωre− Vx
Ωre
. (6)
If the side slip angle is non-zero, lateral slip or side slip is given which becomes, with the used coordinate system and the definition given in [3], as follows
tan(α) = Vy
Vx
. (7)
2.1.2 Tire forces and moments
The horizontal forces, Fx and Fy, and the aligning moment, Mz, are in steady state motion func- tions of the previously introduced slip quantities. Thus slip needs to be present in order for forces
and moments to exist. The wheel load and camber angle has also an influence of the forces and moments.
A turning wheel, under side slip forces causes the treads belonging to the tire to deflect both longitudinally and laterally [4]. This causes an lateral tangential stress distribution and the resul- tant lateral force gets displaced a distance behind the center line [4]. This distance is termed as the pneumatic trail, which gives rise to the aligning moment, Mz. The aligning moment tends to turn the tire towards the direction of the tires velocity vector, V, thus it aligns the wheel with the velocity vector hence its name.
Figure 4 presents typical force curves for pure longitudinal and lateral slip, as well as combined slip. Pure slip occurs when only one slip quantity is present.
Figure 4: Graphs showing typical force characteristics for pure and combined slip [3].
As can be observed in Figure 4, the relationship between the slip and forces is of non-linear nature.
If the slip values are low, which is the case for normal driving conditions [9], a linear relationship between the slip values and forces are valid.
The constant of proportionality for the longitudinal force, Fx, and longitudinal slip, κ, is termed as the longitudinal slip stiffness, CF κ. Similarly the constant of proportionality for the lateral force is termed as the lateral slip stiffness or the cornering stiffness, CF α [3]. Thus for low slip values, following relationships are valid
Fx= CF κκ, (8)
Fy= CF αα. (9)
As further can be observed from the figure, when the situation is such that combined slip is occurring, the forces reduce in magnitude. This behavior is due to that the total horizontal force, F , cannot exceed the maximum friction force, µFz, which is a function of normal force and current surface friction [3]. By the assumption that the tire has isotropic adhesion properties, the maximum allowable horizontal force is given from [10] as follows
Fx2+ Fy2= (µpmg)2. (10)
where µp is the peak coefficient of friction, m the mass and g the gravitational constant. Thus the horizontal force is constrained to be within what’s termed as a friction circle. However this is a rough simplification, and in the reality the tire has anisotropic adhesion properties. Tires which posses the non-linear nature presented in the figure, may be described by a set of empirical formulas called the Magic Formula which are addressed in more detail in section 2.2.1.
2.2 Tire model: Tire models
As was discussed in previous section, the force characteristics of the tire are non-linear. In the most general case, the tire can be modeled as a non-linear system which is dependent on multiple inputs such as its slip quantities, normal load, camber angle etc. The outputs of this general tire model are the tire forces and moments. There exist many different models for tires which couples the previous mentioned quantities. One of the most established tire model is the Magic Formula, which is a semi-empirical tire model and is described in section 2.2.1.
Even though the Magic Formula serves as a satisfactory tire model, it has its limitations to quasi steady state conditions, [11]. Thus when the use case is of transient nature and conditions where the surface is non-flat with short wavelength, the Magic Formula is not applicable. An improved model to the Magic Formula is the MFSWIFT model, which handles transients and is further discussed in section 2.2.2. The brush model is a theoretically derived tire model, and relatively simple in its formulation. The model is described in section 2.2.3. Section 2.2.4 introduces a simple transient tire model, which uses the Magic Formula to calculate the forces.
2.2.1 The Magic Formula
As discussed in previous section, the slip and force relationship for a tire is highly non-linear.
A widely used and well established tire model is the Magic Formula (MF). The MF has been published in many different versions, however the MF used for this thesis work and described in the following section is given from [3], which introduces the general appearance of the MF for the case of pure slip as
y =D sin(C arctan(Bx − E(Bx − arctan(Bx))) (11)
Y (X) =y(x) + Sv (12)
x =X + SH (13)
In the above formula, Y represents the output variable which could either be Fx, Fy or Mz. As input the formula takes the corresponding steady state slip quantitiy, which either is κ or α.
The coefficients in the formula represents a certain characteristic in the curve and some of these coefficients are presented in Figure 5. The coefficients are dependent on normal load, Fz, and inflation pressure, pi. The non-dimensional parameters, dfzand dpi, are introduced in [3], and are as follows
dfz= Fz− Fzo
Fz , (14)
dpi= pi− pio pi
. (15)
The parameters Fzo and pioare the nominal load and the nominal inflation pressure.
Figure 5: The figure illustrates some of the implications that some of the coefficients in the MF has [5].
To explain the relationship between Fz, pi, and the MF coefficients, the longitudinal version of the MF given by [3] will be presented.
The product BCD corresponds to the slope of the linear part of the function. In [3] it’s given a dependency on both load and inflation pressure, and expressed as follows
Kxκ= Fz(PKx1+ PKx2dfz2) exp(PKx3dfz)(1 + ppx1dpi+ ppx1dp2i). (16) The parameter C in Eq (11), controls the range of the sine function which has the implication of controlling the shape of the curve. It’s therefore named as the shape factor, and it can analytically be calculated [3] as
C = 1 ± (1 −2
πarcin(ya
D)). (17)
In the above equation, yais the force which the curve approaches asymptotically as the slip grows.
The coefficient D is the peak value, and controls the peak tractive/braking force of the curve. The peak force is limited by the friction coefficient, µx, and dependent on the normal load. In [3], this is expressed as
Dx= µxFz. (18)
where the variable µxfurther is given a load and inflation pressure dependency as follows
µx= (PDx1+ PDx2dfz)(1 + PP x3dpi+ PP x4dp2i). (19) The factor E is termed as the slope factor, which can be calculated by using the information of the location of the peak force, κm. In [3], a analytical expression for the slope factor is given as
E = Bxκm− tan(2Cπ
x)
Bκm− arctan(Bκm), if C > 1 (20) and its dependency of normal force is in [3] expressed as follows
Ex= (PEx1+ PEx2dfz+ PEx3dfz2)(1 − PEx4sign(κx)). (21) The coefficient Bx is the only free coefficient to determine the initial slope, thus its termed as the stiffness factor and calculated as
Bx= Kxκ/(CxDx). (22)
For some cases, the curve passes through the origin, but this is not always the case since effects such as conicity and ply-steer, which are effects connected to non-symmetry in the tire construction, offset the curve [3]. The rolling-friction also has this effect. To allow the formula to describe this behavior the coefficients SHx and SV xare introduced, which offset the curve horizontally and vertically. They are in [3] expressed as
SHx = (PHx1+ PHx2dfz), (23)
Combined slip
For the case of combined slip (α 6= 0 and κ 6= 0), a weighting factor is introduced [3]. It’s multiplied with the pure slip force to give the effect of a combined slip condition. In [3], following weighting function is presented
G = D cos(C arctan(Bx)). (25)
The coefficients in the function have similar interpretations as the previously introduced. In [3]
the coefficients are presented as follows. D is the peak value, C determines the height and B the sharpness of the function. In the ideal case, the behavior of the function is such that the function equals unity in the case of pure slip which gradually decreases as the situation turns into combined slip [3]. The function therefore takes the form of a hill. Figure 6 and 7 illustrates the weighting function and the effects it has on the lateral and longitudinal force when combined slip is present.
Figure 6: Graphs the weighting function (left graph) and the influence is has on the lateral force (right graph) [3].
Figure 7: Graphs the weighting function (left graph) and the influence is has on the longitudinal force (right graph) [3].
2.2.2 MFSWIFT
The MF is accurate and robust, but is limited to describe the forces under quasi steady state conditions. To describe the dynamics in more detail, the MFSWIFT model could instead be used, which is a fusion of the MF together with the SWIFT model. A conceptual overview of the model is presented in Figure 8.
Figure 8: Illustrates the concept of the MFSWIFT model [12].
The model consists of four important elements [13], which are the MF (described in previous section), a rigid ring to approximate the dynamics of the tire belt, the contact patch slip model, and finally a enveloping model to handle obstacles and arbitrary 3 D roads.
Rigid ring dynamics
To properly describe the dynamics of the wheel, a rigid ring is included into the model. The purpose of this is to include the inertia of belt as well and the tire sidewalls [12]. By including the inertia of the belt, the wheel is split into two parts, rim and a tire belt, which are interconnected with each other by a set of dampers and springs in all directions. The belt is further connected to the contact patch by a residual spring. By dividing the wheel into two parts, the model becomes fairly accurate in the frequency range (60-100Hz [3]), due to that the first mode shapes are rigid vibration modes, which makes the rigid ring approximation valid.
The contact patch slip model
A certain distance is needed for the tire to be travelled before the steady state forces are reached.
The distance needed to be travelled is termed as the relaxation length and can be interpreted as a spatial time constant. From experiments [13] it has been observed that the relaxation length is dependent on the load, magnitude of slip and inflation pressure [14]. An example of this dependency on the side slip is shown in Figure 9, where a step response in side slip is introduced. It’s seen that with increasing side slip the tire responds quicker, thus the relaxation length decreases with increasing side slip.
Figure 9: Illustrates how the magnitude of side slip affects the relaxation length [13].
To describe these tire transients, a contact-patch model which accounts for the transient behavior
analytical model that describes the non steady state response to slip variations is approximated by a set of first-order differential equations, where the transient slip values are given as outputs from these sets of equations. To calculate the transient slip forces, the transient slips are used as inputs in the MF. The forces generated are then applied to the contact patch which further propagates to the wheel axle through the set of dampers and springs, which connects the different parts [13].
Enveloping model
The previous discussion has been limited to a flat and smooth plane road or a road where the wavelength of the vertical road height is sufficiently large. In these situations the tire is excited by axle motion, braking or steering. However the tire can also be excited via the road, examples are road unevenness with short wavelengths (two to three time the contact length [15]) or when traversing obstacles. This gives rise to non-linear behavior of the horizontal and vertical force as well as changes of the effective rolling radius re [13].
To incorporate this excitiation effect into the model, the actual road is filtered by using a set of elliptical cams which are discretized [12] and used to sense the modeled road, this is illustrated in Figure 10.
Figure 10: Presents the concept of elliptical cams [14] . 2.2.3 The brush model
The brush model is a physical tire model, which approximates the tire as a rigid ring where brushes extend radially outward from its circumference [6]. The bristles are given a linear stiffness and are supposed to model the elasticity of the carcass, belt and treads [3]. A schematic layout of the brush model is presented in Figure 11.
Figure 11: Illustration of the brush model [6].
If the tire rolls freely, without any presence of slip, the bristles remain vertical when entering the contact region. No horizontal deflection is given, why no forces are introduced. If however slip is present, the brushes deform horizontally in the contact area which generates the tire forces and moments [3].
In the contact area, two regions are defined, which are the adhesive and sliding region [3]. Within the adhesive region the bristles follow the direction dictated by the velocity vector, which by def- inition is the side slip angle α. For the case of pure lateral slip, the tip of the bristle enters the
contact region with zero lateral deflection which then grows linearly as the tip progresses through the contact area. The maximum possible deflection is however limited by the present friction co- efficient, the vertical force distribution and the stiffness of the element [3]. A transition from the adhesive region to the sliding region occurs when the limit is reached, and the deflection instead follows a pattern dictated by the previous mentioned quantities. This is illustrated in Figure 12 for the case of a parabolic pressure distribution. If the situation is such that pure longitudinal slip occurs, the bristles move rearwards in the case of braking and forward in the case of driving [6].
Figure 12: The left illustration shows how the bristles deflects as side slip increases. The plot to the right plots the corresponding lateral forces [3].
Pure longitudinal slip
The brush model in a state of braking and under pure longitudinal slip is illustrated in Figure 13.
Figure 13: The brush model for the case of pure longitudinal slip [3].
In the model, the base points of the bristles move rearwards with the velocity Vr termed as the linear speed of rolling, which is as follows
Vr= Ωre. (26)
The elements tip points adheres to the ground which results in a longitudinal deflection if slip is present. In [3] this is expressed as
u = −(ac− x) Vsx
Vx− Vsx
. (27)
In [3], the theoretical longitudinal slip, σx, is introduced and defined as σx= −Vsx
Vr
= κ
1 + κ. (28)
By the introduction of the theoretical slip, the longitudinal deflection may instead be expressed in terms of this quantity, which by [3] becomes
u = (ac− x)σx. (29)
A composite function θxis in [3] introduced θx= 2
3 cpxa2c
µFz
(30) and total sliding occurs at
σslx= ± 1
θx. (31)
If a parabolic force distribution of the normal load is assumed in the contact area, then following expression for the longitudinal force is given [3]
Fx=
3µFzθxσx(1 − |θxσx| +1
3(θxσx)2) , if |σx| ≤ σslx
µFzsgn(κ) , if |σx| ≥ σslx.
(32)
2.2.4 Single contact point model
The single contact point model (SCP) [3] is a transient tire model, which works together in com- bination with the steady state MF to obtain the tire forces. It’s limited to low levels of slip [11].
The model is illustrated in Figure 14.
Figure 14: Illustration of the single contact point model [3].
From the description given in [3], two different points are used, the slip point S which is attached to the rim and the contact point S0. These points move with different velocities, and the difference between them give a longitudinal deflection, u, and lateral deflection, v. The differential equation for the longitudinal case reads according to [3] as
du dt +1
σ|Vx|u = −|Vx|κ = −Vsx (33)
where σ is the longitudinal relaxation length. Next, the transient slip quantity is calculated [3] as κ0 =u
σ. (34)
This transient slip quantity is then used as input in the MF in order to get the forces acting on the tire
Fx= M F (κ0, Fz) (35)
where M F is the Magic Formula for the pure longitudinal slip case.
2.3 Tire model: Algorithms for parameter estimation
The MF is a widely used and established tire model, but it possesses a certain level of complexity due to its non-linear relationship and vast number of parameters that it contains. These param- eters are usually determined by least square fitting techniques [7]. The techniques uses what is termed as a cost function, which defines how good of a fit some parameters give to observed data.
The problem then becomes to find the parameters which gives the best fit.
In [7] various algorithms were applied to estimate the parameters from the MF. The used al- gorithms, included both gradient based and derivative-free methods as well as unconstrained and constrained methods. The algorithms were applied to lateral force, longitudinal force and self- aligning torque data, all under pure slip conditions. The data that were used was also from
different sources. The three methods that gave best fits were the trust region reflective, differential evolution and bounded cuckoo search. A important conclusion from [7] was that even though an algorithm converges to a certain set of parameters, the calculated parameters can differ between different algorithms. A recommendation from [7] was to compare the results of a number of algo- rithms before choosing the parameters. If the estimated parameters are all the same, then clearly this is a good indication.
The differential evolution algorithm showed good results in [7], why the algorithm was deemed as being suitable for the thesis work. The following section gives a brief description of how the algorithm works.
2.3.1 Differential evolution
Differential evolution is based on evolutionary optimisation techniques [7], which means that it mimics the way evolution works. The algorithm works by defining a initial population, Gi. Each inhabitant is a vector and its elements contains the parameters that are suppose to be estimated.
These are in the beginning chosen randomly [16]. A new generation is created by a process termed as mutation [16]. Practically this is performed by combining randomly chosen inhabitants in the current generation, to create a new mutated vector. As a next step, the mutated vectors parameters are mixed with a chosen inhabitant in the current generation, Gi, to give diversity in the solutions [16]. This process is termed as crossover. By performing these steps, a new generation is given, Gi+1, which contains new potentially better parameters. The process of deciding if a a inhabitant in the new generation, Gi+1, should replace a corresponding inhabitant in the current generation, Gi, is termed as the selection process [16]. The previously mentioned steps are then iterated until a solution is found.
2.4 Force distribution: Normal force
The theory regarding normal force distribution and estimation is quite extensive for the automo- tive industry. There are some papers regarding articulated machines [17, 18] but when looking for articles on modelling articulated machines most articules focuses on trucks. The articulation system for wheel loaders differ quite a bit from trucks as the wheel loaders articulation joint is unable to turn by itself as it is held stiff and controlled by two hydraulic cylinders. However when including tractors and general articulated machines as a base for modelling wheel loaders there are some articlesthat are useful [19, 20, 8, 21]. These paper however do not describe the load distribution in a way that is useful for this thesis. The limit on any model created is to only use available sensors which can measure: steering angle, acceleration of the machine, rotational wheel speed, motor torque, change in the machine articulation angle, angle of the lift arm and lift arm hydraulic cylinder pressure caused by the load. Anything else is considered unavailable and any model requiring more inputs than these, unusable. As such it was determined to build a model from the ground up, using simple physics in order to have a better understanding of the model and be able to improve it as time went on. When deciding how to model the system, inspiration was collected from several sources, including models for cars [22, 8, 20]. Due to the available sensor a local coordinate system was chosen, this means the machine only know its position with regards to its own origin point placed on the machine, but has no idea about its position relative the global world. General vehicle dynamics effects were also investigated [23, 3, 24] when building this model.
With the introduction of powerful electric motors the wheel loaders are able to perform in ways not previously possible. This creates situations that were previously safe which can now, with for example faster acceleration, lead to roll over. For inexperienced drivers it becomes extra important to prevent these situations from arrising. When it comes to roll-over prediction it has been pro- posed that looking at the difference between left and right side wheel load can be used to estimate when roll over is likely to happen [8, 18] so this will be the approach for this thesis as well.
2.4.1 Defining the basis for the normal force distribution model
How the center of gravity of a machine affects the normal force on the tires is the basis for the model to be used. When calculating the distribution, the origin is placed at the articulation joint.
The x-direction is defined positive in the front-facing direction of the machine, the longitudinal direction. The y-direction is positive to the right of the machine, in the lateral direction and the z-direction is positive upwards, in the vertical direction. Figure 15 shows a simple illustration of a wheel loader from above with the coordinate system set at the articulation point.
Figure 15: Wheel loader as seen from above facing ”up-ward”.
2.4.2 Finding the center of gravity
The wheel loader is able to turn at the middle due to its articulation joint. In order to estimate the normal forces on each tire it is important that the static center of gravity is accurately calculated.
The static center of gravity refers to when the machine is in a steady state and accounts only for the articulation angle δ, the lift arm angle φ and the bucket load. The machine is divided into separate parts: The rear, the front, the lift arm and bucket. The lift arm and bucket will collectively be named simply as bucket from here on out. The center of gravity is calculated independently in the x, y and z directions by taking how far away the rear, front and bucket center of gravity is from the articulation point and weighting the values with the parts different masses. Figure 16 shows the distance in the x-direction for a machine that is turning.
Figure 16: Distance between masses and origin in the x-direction.
The machine center of gravity is calculated in the following way
Xcog =
mf rontxf rontcos(δ) + mbucketLoadxbucketLoadcos(δ) + mbackxback 1 mmachine
. (36) Xcog is the distance in the x-direction from the articulation point to the center of gravity. m represents mass and x the longitudinal distance for the given part. The weight of the bucket with
arm-linkage and the payload has been combined, while accounting for the fact that the change in arm position will happen some distance from the origin. This then becomes
xbucketLoad= 1 mbucketLoad
(mbucket(xf ront+ (xbucket− xf ront) cos(φ) +mload(xf ront+ (xload− xf ront) cos(φ)).
(37)
mbucketLoad is the combined weight of the bucket, arm-linkage and load
mbucketLoad= mbucket+ mload (38)
and mmachineis the total mass and is the sum of the individual masses
mmachine= mback+ mf ront+ mbucketLoad. (39) In the y-direction it becomes
Ycog= (mf rontxf rontsin(δ) + mbucketLoadxbucketLoadsin(δ) + mbackxback) 1 mmachine
. (40) Since the contact is along the ground the vertical distance will be of importance for the moment equilibrium equations and is
Zcog= (mf rontzf ront+ mbucketLoadzbucketLoad+ mbackzback) 1 mmachine
(41) and
zbucketLoad=
mbucket(zbucket+ (xbucket− xf ront) sin(φ)+
mload(zload+ (xbucket− xf ront) sin(φ) 1 mbucketLoad.
(42)
The distribution of the load on the wheels is assumed to be inversely proportional to the distance from the center of gravity. The first step is then to find the position of the wheels , wp, with regard to the origin. These are calculated as vectors for the x- and y-coordinates to make the equations more clean. The wheel loader is symmetric for the wheels along the x-axis when straight, so the following equations hold true:
wpf ront lef t= xf rontAxlecos(δ) − yf rontAxlesin(δ)
−yf rontAxlecos(δ) − xf rontAxlesin(δ)
, (43)
wpf ront right=xf rontAxlecos(δ) + yf rontAxlesin(δ) yf rontAxlecos(δ) − xf rontAxlesin(δ)
, (44)
wprear lef t= xrearAxle
−yrearAxle
, (45)
and
wprear right=xrearAxle yrearAxle
(46) where the longitudinal distances are between the articulation point and the middle of the axles and the lateral distances are to the end of the axle, at a point in the middle of the tire.
2.4.3 Shifts in center of gravity
To account for shift in the center of gravity due to change in speed of the machine, road inclination or rocking of the machine, equations of equlibrium will be used. As the normal force estimator will be an input to the tire model it can not use the longitudinal or lateral forces, since those are calculated in the tire model. Instead the acceleration will be used since that can be measured or calculated from avaliable sensors. For the change in angle with regards to inclination of the road or rocking of the machine, this data is avalaiable from gyroscopes installed on the machine. To calculate the shifts in center of gravity, the theory used for cars will be used as the base and then later adopted for the center of gravity model. All equlibrium equations are derived from [4].
Steady-state longitudinal case
A stationary machine is a simple system but will be used to define the system used for weight distribution calculations. The system is shown in Figure 17.
Figure 17: Stationary machine with the front to the right.
The longitudinal distance between wheel position and center of gravity is represented by a and index f refers to the front of the machine while b refers to the back of the machine. The forces on the machine when stationary can be described by two equilibrium equations
XFz= 0, XMy = 0.
(47)
In this case it is assumed that the weight is roughly in the middle of the left and right side. The sum of the normal forces will equal the weight of the machine
XFz= 2Fzf+ 2Fzb− mg = 0. (48)
The moment around the y-axis comes from the vertical forces acting on the tires
XMy= 2Fzfaf− 2Fzbab= 0. (49)
Rewriting the equation to find the relation between the vertical forces yields Fzb= Fzfaf
ab. (50)
Substituting into Eq (48) gives
2Fzf+ 2Fzfaf
ab − mg = 0. (51)
Simplifying and solving for Fzf produces
Fzf =1
2mg ab
af+ ab
. (52)
The force on the back tire is then calculated from Eq (50).
Accelerating in a straight line
The next step is to look at the how the vertical forces shifts between the tires when the machine is accelerating. The forces are shown in Figure 18.
Figure 18: Accelerating machine with the front to the right.
This time the equilibrium equations become
XFx= m ˙v, XFz= 0, XMy= 0.
(53)
The longitudinal forces will be equal to the acceleration times the mass of the machine
XFx= 2Fxf+ 2Fxr= m ˙v (54)
where v is the velocity of the machine and ˙v is its time derivative. The sum of the normal forces will still be the total weight of the machine
XFz= 2Fzf+ 2Fzb− mg = 0. (55)
This time the moment around the y-axis is caused both by the vertical forces at the tires, but also by the longitudinal forces
XMy= −2Fzfaf+ 2Fzbab− (2Fxf+ 2Fxr)Zcog= 0. (56) Substituting Eq (54) into Eq (56) gives
XMy= 2Fzfaf− 2Fzbab+ m ˙vZcog = 0. (57) Solving for Fzb
Fzb =m ˙vZcog+ 2Fzfaf 2ab
(58) and substituting into Eq (55) gives the solution for Fzf
2Fzf+ 2m ˙vZcog+ 2Fzfaf 2ab
− mg = 0, Fzf =1
2
mg ab af+ ab
− m ˙v Zcog af+ ab
.
(59)
The back tire force is calculated from Eq (58).
Stationary machine on inclined road
When looking at the case of a machine on an inclined road, the stationary case is shown in Figure 19.
Figure 19: Inclined machine with the front to the right.
The angle θ shown in the figure is the pitch angle. The equilibrium looks similar to the accelerating case, except the sum of the longitudinal forces are zero
XFx= 0, XFz= 0, XMy = 0.
(60)
The longitudinal forces must resist the part of the weight which is trying to move the machine backwards
XFx= 2Fxf+ 2Fxr− mg sin(θ) = 0. (61) The vertical forces only takes a portion of the total weight, the part that is pushing against the ground
XFz= 2Fzf+ 2Fzb− mg cos(θ) = 0. (62) The moment around the y-axis is by itself independent of the angle of the incline. Any variance in the moment comes from the fact that the effecting forces change with the angle
XMy= 2Fzfaf− 2Fz2ab+ (2Fxf+ 2Fxr)Zcog = 0. (63) Inserting Eq (61) into (63) and simplifying in the same way as before gives
Fzf =1 2mg
af
af + ab cos(θ) − Zcog
af+ absin(θ)
, Fzb =1
2mg
ab af + ab
cos(θ) + Zcog af+ ab
sin(θ)
.
(64)
Stationary machine on banked road
A shift of the weight in the lateral direction can happen for example when the machine is on a banking road as shown in Figure 20.
Figure 20: Machine on a banking road, seen from behind.
Index r refers to the right side and l to the left while ψ is known as the roll angle and b is the lateral distance between wheel and center of gravity. For the equilibrium equations the moment will be calculated around the y-axis instead of the x-axis
XFx= 0, XFz= 0, XMy = 0.
(65)
The lateral forces counteract the lateral portion of the weight which can be seen in
XFy= 2Fyr+ 2Fyl− mgsin(ψ) = 0. (66)
The vertical forces correspond the the part of the weight that is directed into the inclined ground as shown in
XFz= 2Fzr+ 2Fzl− mgcos(ψ) = 0. (67)
The moment around the y-axis is positive in the clock-wise direction as follows
XMx= 2Fzlbl− 2Fzrbr− (2Fyr+ 2Fyl)Zcog= 0. (68) The forces on each tire is calculated in the same way as prevision examples resulting in
Fzr= 1
2mg bl br+ bl
cos(ψ) −1
2mg Zcog br+ bl
sin(β) Fzl= 1
2mg br br+ bl
cos(ψ) +1
2mg Zcog br+ bl
sin(ψ).
(69)
2.4.4 Normal force on each tire
Adapting the equation above for an articulated wheel loader is possible when all distances between the wheels and the center of gravity are determined. To reduce the size of the equations and make them easier to follow each of the four tires will be given a separate numerical index and the mass will be designated mt. The numeration convention is as follow
1. Front left tire.
2. Front right tire.
3. Back left tire.
4. Back right tire.
Thich is illustrated in Figure 21
Figure 21: Labelling of wheels on a wheel loader.
To find the forces on each tire, five equilibrium equations can be used to evaluate the effect off the normal forces
XFz= 0, XFx= mt˙vx, XFy= mt˙vy, XMy= 0, XMx= 0.
(70)
But before those equations are defined it is necessary to look at the rear axle. As there is no damping in a wheel loader the rear axle can rotate independently of the rear body, which is where the cabin and driver are located. This changes the way the load is distributed. Figure 22 shows a detailed view of a wheel loaders rear axle.
Figure 22: Rear axle that can rotate independent of the body [25].
The mounting points are attached to the wheel loader body and join together to create a single contact point for the rear. This means that in order to get a accurate machine model, the wheel loader has to be modeled as a tricycle. This set up is shown in Figure 23.
Figure 23: Tricycle configuration of load distribution.
Using the tricycle model where mt is the mass of that portion of the machine, the sum of the vertical forces becomes
XFz= Fz1+ Fz2+ FzR− mtg cos(ψ) cos(θ) = 0 (71) where index R refers to the forces on the rear contact point. The longitudinal forces are
XFx= Fx1+ Fx2+ FxR− mtg sin(θ) = mt˙vx (72) and the lateral forces are
XFy= Fy1+ Fy2+ FyR− mtg sin(ψ) = mt˙vy. (73) Before calculating the moment around the x- and y-axis it is important to look at where the rear forces act. The tire forces act on the ground at the contact points between the wheels and ground, but the rear forces acts where the support attaches to the rear body, as can be seen in Figure 22.
The configuration instead becomes the one seen in Figure 24.
Figure 24: Contact points for the various forces for the tricycle model.
With this in mind the moment around the y-axle is
XMy= Fz1a1+ Fz2a2− FzRaR+ Zcog(Fx1+ Fx2) + (Zcog− ZR)FxR= 0 (74)
