Longitudinal and Lateral brush tyre model using Matlab

Tyre Modeling

A tyre model is a mathematical or an empirical representation of the mechanism through which the tyre generates forces and moments in response to inputs such as slip angle and camber angle

Brush Tyre Model

The elastic band model or the brush tyre model helps to explain the relationship between Longitudinal force (`F_x`) and Longitudinal slip (slip ratio) using fore-aft deflection of the brushes. This is later extended to explain the relationship between Lateral force (`F_y`) and slip angles by considering the lateral deflection of the brushes. In the combined slip model, both the longitudinal and lateral forces and slips are studied together along with the effects they have on each other. The Brush tyre model is based on the principles of Julien's Theory which makes the following assumptions about the tyre:

• Tyre threads are considered to be elastic bands or brushes (springs) instead of the usual consideration as blocks.
• The contact patch is considered rectangular in spite of the fact that due to rolling resistance, the pressure distribution along the contact patch is irregular giving it an asymmetrical shape.
• The pressure distribution in the contact patch is uniform
• The contact patch is divided into adhesive and sliding regions.

Slip Ratio and Longitudinal Force:

In order to formulate the longitudinal brush tyre model, it is necessary to understand the generation of slipping due to the application of longitudinal force. Longitudinal slip is generated due to braking and acceleration forces. The mechanism of slip generation is the same for both the cases where a difference between the velocity of the tyre at contact patch and the velocity centre progressively gives rise to slip.

Braking Force

When brakes are applied, a negative torque known as braking torque opposes the relative motion of the wheel. This application of brake torque (`T_b`) leads to a difference between the velocity at the contact parch and the velocity at the wheel centre. This is because, at the contact patch, the tyre experiences an angular velocity due to braking (`omega_b`), which is less than the actual velocity (`omega`). The difference in velocities slides the thread material relative to the road and gives rise to slip. The braking force pushes the contact patch rearward causing tension at one end compression at the other. As the longitudinal slip increases the shear stress builds up until the point where it cannot exceeds the limit of friction and starts to slide. This sliding happens at the rear of the contact patch. If the brake force is maintained, the angular velocity drops from its free rolling value eventually to zero when the wheel is locked.

Longitudinal shear force distribution

1. Free rolling 2. Braking 3. Net result

Slip ratio `S_x` is due to the difference in velocities of the free-rolling wheel (`omega_o`) and braked wheel (`omega_b`)

`S_x = (omega_o-omega_b)/(omega_o)`

for free rolling forward velocity of the wheel

`v = omega_o * R_e`

`omega_o = v/R_e`

`S_x = (v - (omega_b*R_e))/v`

when the wheel is locked then `omega_b` = 0 and `S_x = 1` then it causes 100% longitudinal slip.

For a given vertical load (`F_z`), the braking force increases linearly with slip angle up to a point (typically 15 to 30 %) beyond which it saturates. ABS (Anti-lock braking system) systems cycle the brake pressure to maintain the slip ratio where maximum braking occurs in order to maximize braking effort and also to maintain a rolling wheel. The peak brake force depends strongly on the tye road interface(coefficient of friction `mu`)

Driving ( Acceleration) force:

During acceleration, the threads get compressed on the front end and extend (tension) at the rear end because the driving force pushes the contact patch rearward. The resultant vertical force (`F_z`) moves closer to the tyre centre. The transmission imparts a driving torque (`T_d`) at the rotating wheel. This leads to a difference in velocities between the tyre contact patch and wheel centre. As a result, the thread elements begin to slide relative to the road. As the thread elements enter and move forward, the longitudinal shear stress builds up in a linear manner.

1. Free rolling 2. Driving 3. Net result

`S_x = (omega_d - omega_o)/omega_d`

`v = omega_o * R_e`

S = `(omega_d *R_e - v)/(omega_d *R_e)`

if v = 0 ; vehicle is stationary and wheel is spinning then S = 1 or 100%

Traction control (TC) systems are used to maintain the slip ratio 15 to 30 % to get the formal propulsive force.

Brush Tire Model - Equations

Longitudinal Brush tyre Model:

This model has been formulated for acceleration slip and it can be modified to include braking slip too. However, the model does not include the sliding zone as it is highly non-linear and can become difficult to formulate. However, the effect of the tyre during translational zone can be extended to show the sliding zone of the tyre. The formulation is done for three different regions of the contact patch:

• Adhesion region where critical slip length is greater than contact patch length indicating that the tyre is in the linear zone.
• The critical region where critical slip length is equal to the length of the contact patch. This is the end of the adhesive region where the relationship between normal force and the maximum longitudinal force is derived.
• The translational region where the tyre in the adhesive region and partly in the sliding region. Beyond the peak value of traction force, the tyre goes into the sliding zone.

`l_t` `to` length of contact patch

`omega_r` `to` angular velocity of the tyre

Strain in the contact patch = `Delta_l/l = ((R*omega_r* Delta_T) - (v*Delta_T))/(R*omega_r* Delta_T)` = slip

brush deformation = `Delta_l` = Strain `(or slip) * l`

`Delta_l = R*omega*Delta_T*i`

e = `Delta_l `

x = `R*omega*Delta_T`

i = slip ratio

e = x*i

Initial deformation `e_o` = `lamda *i`

Total deformation = `(lamda+x)*i`

Force per unit length is directly proportional to Total deformation.

`d/dx(F_x) = k_t*(lamda+x)*i`

`k_t` `to` tangential stiffness (k-N/m^2)

For adhesion, `d/dx(F_x) <= (mu_p*F_z)/l_t to mu_p*P*b`

`mu_p` to peak friction coefficient

`P` `to` Pressure

`b` `to`width of the contact patch

For Critical condition , `d/dx(F_x) = (mu_p*F_z)/l_t = k_t*(lamda+x)*i`

`x = (mu_p*F_z)/(l_t *k_t*i) - lamda`

`x to` critical length `(l_c)` `to` distance along contact patch

if `x < (mu_p*F_z)/(l_t *k_t*i) - lamda` `to` Adhesion

`x > (mu_p*F_z)/(l_t *k_t*i) - lamda` `to` Sliding

Critical slip `i_c`

`k_t*(lamda+x)*i` = `(mu_p*F_z)/l_t`

`i_c = (mu_p*F_z)/(l_t*k_t*(lamda+x))`

1) Adhesion region

slip < critical slip

`l_c > l_t`

`int_0^(l_t) (dF_x) = int_0^(l_t) k_t*(lamda+x)*i *dx`

= `k_t*lambda*i*x+ k_t*i*x^2/2`(limits 0 to `l_t`)

`implies F_x = k_t*l_t*i*[lambda +l_t/2]`

special case: if `lambda` = 0

`F_x = 1/2*k_t*l_t^2*i`

`C_i = 1/2*k_t*l_t^2 to` longitudinal stiffness

`F_x = C_i *i`

2) Critical condition

`F_(xc)` to end of linear zone

`F_(xc) = k_t*l_t*[[2*lambda + l_t]/2] *mu_p*F_z/(k_t*l_t*(lambda+l_t)) `

`F_(xc) = mu_p*F_z/((lambda+l_t))* (2*lambda+l_t)/2 `

if `lambda = 0`

`F_(xc) = mu_p *F_z/2`

3) Sliding

`F_(xa)` (adhesion) = `int_0^(l_c) k_t*(lamda+x)*i *dx implies k_t*l_t*i*[lambda +l_c/2]`

`F_(xs)` (sliding) = `int_(l_c)^(l_t)(mu_pF_z)/(l_t)dx implies (mu_pF_z)/(l_t)[1-l_c/l_t] `

`F_(x)`total = `F_(xa) + F_(xs) = k_t*l_t*i*[lambda +l_c/2] + (mu_pF_z)/(l_t)[1-l_c/l_t] `

Special case `lambda` = 0

`l_c = (mu_p*F_z)/(l_t *k_t*i) - lamda`

`F_(x)`total = `(k_tl_c^2i)/2+mu_p*F_z - mu_p*F_z*l_c/l_t`

`F_x = mu_pF_z - (mu_pF_z)^2/(4C_ii)` where `C_i = (k_t*l_t^2)/(2)`

Three parameters are required for Longitudinal brush tyre model

`F_z , mu_p, C_i`

Braking

`i_s` `to` skid ratio

Slip `to` forward acceleration

Skid `to` braking

`i_s = (v-omega *r)/v`

`i_s/(1-i_s)`

`((v-omega *r)/v)/(1-(v-omega *r)/v)` `implies` `(v-omega*r)/(omega*r)`

The equations formulated above can be modelled in Matlab/Octave/Python to see the effect of the longitudinal force in slip generation. The modelling was done to study the effect of the peak friction co-efficient on Traction Force generation and critical length. The code would take user-specified inputs for tyre parameters like normal force, tyre tangential stiffness, peak friction co-efficient `(mu_p)`, contact patch length and initial contact patch deformation.

From the below plots, it can be seen that the critical length is reducing with an increase in longitudinal slip and eventually it reaches zero at a point. At this point the traction force is maximum. This is valid only when the contact patch begins to slide because sliding friction is not included in the above equation. Also, as the coefficient of friction reduces its peak traction force also reduces and the critical length reaches to zero faster as compared to the former case.

Matlab/Octave code for Longitudinal Brush Tyre Model

clear all
close all
clc

% To_plot_the_relationship_between_traction_force_slip_ratio
% To_plot_the_relationship_between_critical_length_and_slip_ratio
%input data

Fz = 4000; % Normal_force (N)
mu = 0.8;  % mu_p = peak_road_tire_friction_co-efficient
lt = 0.2;  % contact_patch_length (m)
kt = 3E6;  % tire_tangential_stiffness (N/m^2)
lambda = 0.04; % Initial_deformation

%compute_critical_slip
Sx_c = mu*Fz/(kt*lt*(lt+lambda));

% main loop to compute traction force

for i = 1:101
    Sx(i) = (i-1)*0.01; %Sweep sx from 0 to 1
    
    if Sx(i) < Sx_c %Adhesion_Region
        Fx(i) = kt*lt*(lambda+lt/2)*Sx(i);
        lc(i) = lt;
    else
        lc(i) = mu*Fz/(lt*kt*Sx(i))- lambda;
        Fx(i) = kt*Sx(i)*(lc(i))*(lambda+lc(i)/2)+mu*Fz*(1-lc(i)/lt);
    end
end

subplot(2,1,1);
plot(Sx,Fx,'b','Linewidth',3);
grid on;
set(gca,'fontsize',18);
title('Traction Force vs Slip Ratio');
xlabel('Longitudinal slip Sx');
ylabel('Traction Force Fx');
hold on;


subplot(2,1,2);
plot(Sx,lc,'r','Linewidth',3);
grid on;
set(gca,'fontsize',18);
title('Critical length vs Slip Ratio');
xlabel('Longitudinal slip Sx');
ylabel('Critical length (m)');
hold on;


    
%%% Reduced_Friction_Case_%%%


Fz = 4000; % Normal_force (N)
mu = 0.5;  % mu_p = peak_road_tire_friction_co-efficient
lt = 0.2;  % contact_patch_length (m)
kt = 3E6;  % tire_tangential_stiffness (N/m^2)
lambda = 0.04; % Initial_deformation

%compute_critical_slip
Sx_c = mu*Fz/(kt*lt*(lt+lambda));

% main loop to compute traction force

for i = 1:101
    Sx(i) = (i-1)*0.01; %Sweep sx from 0 to 1
    
    if Sx(i) < Sx_c %Adhesion_region
        Fx(i) = kt*lt*(lambda+lt/2)*Sx(i);
        lc(i) = lt;
    else
        lc(i) = mu*Fz/(lt*kt*Sx(i))- lambda;
        Fx(i) = kt*Sx(i)*(lc(i))*(lambda+lc(i)/2)+mu*Fz*(1-lc(i)/lt);
    end
end

subplot(2,1,1);
plot(Sx,Fx,'m','Linewidth',3);
grid on;
set(gca,'fontsize',18);
title('Traction Force vs Slip Ratio');
xlabel('Longitudinal slip Sx');
ylabel('Traction Force Fx');
legend('mu = 0.8','mu = 0.5');


subplot(2,1,2);
plot(Sx,lc,'k','Linewidth',3);
grid on;
set(gca,'fontsize',18);
title('Critical length vs Slip Ratio');
xlabel('Longitudinal slip Sx');
ylabel('Critical length (m)');
legend('mu = 0.8','mu = 0.5');

Slip Angle and Lateral Force: When a tyre print distorts, it continuously loads and unloads. During loading the energy is higher, whereas when the tyre print unloads it follows a different curve suggesting a hysteresis loss. This energy loss is what gives rise to rolling resistance in tyers where the loaded print has higher energy than that of the unloaded portion. When the cornering force is applied the inertial force (centrifugal force) will resist the cornering effect which in turn causes distortion at the tyre footprint. Friction is what sustains this distortion and leads to the build-up of lateral force/shear stress. Due to this distortion, the wheel travels in a different direction from where it points.

When a rolling pneumatic tyre is subjected to lateral force, an angle known as slip-angle is created between the direction of heading and wheel travel. When the lateral stress build-up caused by the creation of slip angle exceeds the boundary limits, sliding occurs and the lateral stress follows the boundary limit distribution. The area of stress build-up is a measure of the resulting cornering force `(F_y)`. As the shape of the stress distribution is approximately triangular, the lateral force acts through the centroid. This distance between this centroid and the wheel centreline is called the pneumatic trail. This mechanism is not instantaneous and the lateral force generation lags the slip angle by a certain distance and the tyre has to roll called the relaxation length. Sidewind and road disturbances can also cause lateral force.

Lateral Brush Tyre Modelling

Similar to longitudinal brush tyre model, the lateral tyre model also does not include the formulation of the sliding zone as it is highly non-linear. However, the effect of the tyre during translational zone can be extended to show the sliding zone of the tyre. The formulation is also done for three different regions of the contact patch.

• The adhesive region where critical slip length is greater than contact patch length indicating that the tyre is in the linear zone
• The critical region where critical slip is equal to the length of the contact patch. This is the end of the adhesion region where the relationship between normal force and maximum lateral force or friction force is derived.
• The translational region where the tyre is partly in the adhesion region and partly in the sliding region. Beyond the peak value of lateral force, the tyre goes into the sliding zone.

Along with formulating lateral force at each of these zones, the aligning moment given by `M_z` that the tyre generates is also formulated. The Aligning moment is generated because of the shape of stress distribution is approximately triangular leading to the lateral force acting at a point behind the wheel centre. This causes a moment about the wheel centre that tries to straighten the tyre during a turn.

a) Adhesion Region : Lateral force, `F_y = C_(alpha)*alpha`

where `C_(alpha)` is the cornering stiffness, `alpha` is the slip angle

`C_(alpha) = (k_y*l_t^2)/2`; `k_y` is the lateral stiffness of the tyre.

Aligning Moment, `M_z = C_(alpha)*l_t*alpha/6`

b) The critical condition for Peak Adhesion:

Here, critical slip angle, `tan(alpha_c) = mu_p*F_z/2*C_(alpha)`

Critical length, `l_c = (mu_p*F_z)/(2C_(alpha)`

C) Translational Zone (sliding):

Lateral Force, `F_y = mu_p*F_z - (mu_p*F_z)^2/(4*C_(alpha)*tanalpha`

Aligning Moment, `M_z = (mu_pF_z)^2/(12*l_t*k_y*tanalpha)[3-(2mu_pF_z)/(l_t^2*k_y*tanalpha)]`

The equations formulated above can be modelled in Matlab/Octave/Python to see the effect of lateral force in the slip angle generation. The modelling was done to study the effect of the peak friction co-efficient on lateral force and Aligning moment. The code would take inputs for tyre parameters like normal force, tyre lateral stiffness, peak friction coefficient,contact patch length and initial contact patch deformation.

clear all
close all
clc

% Lateral brush tyre model
Fz = 5000; % Normal force(N)
mu = 0.9; % Peak road - tyre friction coefficient
Ca = 50000; % Lateral stiffness (N)
lt = 0.2; % Assumed length of tyre contact patch
Ky = 2*Ca/(lt^2);

% Compute critical slip
alpha_c = mu*Fz/(2*Ca);

% Main loop to compute tyre lateral force

for j = 1:300,
    alpha(j) = (j-1)*0.1*pi/180; % slip angle between 0 and 15 degrees
  if alpha(j) < alpha_c
      Fy(j) = Ca*alpha(j); % Adhesion only
      Mz(j) = Ca*lt*alpha(j)/6;
  else
      Fy(j) = mu*Fz - (mu*Fz)^2/(4*Ca*tan(alpha(j)));
      Mz(j) = (mu*Fz)^2/(12*lt*Ky*tan(alpha(j)))*(3-(2*mu*Fz/(lt^2*Ky*tan(alpha(j))))); % sliding
  end
end

subplot(1,2,1);
plot(alpha*180/pi, Fy);
set(gca,'fontsize',18)
grid on;
xlabel('Slip angle (deg)');
ylabel('Lateral force Fy (N)');
title('Lateral Force vs Slip Angle');


subplot(1,2,2);
plot(alpha*180/pi, Mz);
set(gca,'fontsize',18)
grid on;
xlabel('Slip angle (deg)');
ylabel('Mz (N-m)');
title('Aligning Moment vs Slip Angle');