Pacjeka magic formula for curve fitting and combined slip correction using Matlab

                                    Pacjeka Magic Tyre model

1. Pacjeka Magic Tyre model was developed by Prof.Hans Bastiaan Pacejka and it is widely used to represent empirically and interpolate previously measured tyre force and moment curves. This model is not a predictive model.
2. It is widely used for simulation work and it is one of the most well-established models.
3. The general magic formula has been developed using mathematical functions that relate:  

1.  Lateral force to Slip angle
2. Aligning moment to slip angle
3. Longitudinal force to longitudinal slip ratio           

                                                 The Magic Formula     

                             `Y(X) = D*sin(Carctan(Bx-E(Bx - arctan(Bx))))`

                                            `Y(X) = y(x) + S_v`

                                                `x = X + S_h`

                           Graphical representation of the magic formula

Whereas

 `Y(X)` = Cornering force, braking force or self-aligning moment

`X` = Slip angle or skid ratio

`B` = Stiffness factor

`C` = Shape factor that controls stretching in the x-direction.

`D` = Peak value or Peak factor

`X_m` = Slip angle at peak value 

`E` = Curvature factor that affects transition in the curve and the position `X_m` at which the peak value occurs

`S_v and S_h` = Offsets to account for camber thrust, Conicity, Plysteer or rolling

`y_s` = Steady state value

                                                 `Slope = B*C*D`

                                                       `B = (Slope)/(C*D)`

values of B

 Lateral force - 1.3

 Longitudinal force - 1.65

 Aligning moment - 2.4

                                                   `y_s = D*sin((pi*C)/2)`

                                                        `E` = `(BX_m - tan((pi*C)/2))/(BX_m - atan(BX_m)`

clear all
close all
clc
load tire_1
a =tire_1(:,1); % slip angle in degree
Mz=tire_1(:,2); % Aligning moment in N-m


% plot test data
plot(a, Mz,'o')
set(gca,'fontsize',18);
title('Pacjecka magic formula');
xlabel('slip angle'); 
ylabel('Self-aligning torque Mz')
hold on
Bz=0.2083;Cz=2.17;Dz=60;Ez=-2.4;Shz=0;Svz=0;
for i=1:101
  a1(i)=(i-1)*0.14;
  phi_z=(1-Ez)*(a1(i)+Shz) + Ez/Bz*atan(Bz*(a1(i)+Shz));
  Mz1(i)=Dz*sin(Cz*atan(Bz*phi_z))+Svz;
end
plot(a1,Mz1,'r','linewidth',3), axis([0 20 0 75]);
grid;  
hold off
sum_error=0;
for i=1:size(a,1);
   phi_z=(1-Ez)*(a(i)+Shz)+Ez/Bz*atan(Bz*(a(i)+Shz));
   Mz_h(i)=Dz*sin(Cz*atan(Bz*phi_z))+Svz;
   sum_error=sum_error+(Mz_h(i)-Mz(i))^2;
end
rms_error=(sum_error/size(a,1))^0.5

rms_error =

1.9097

                                                          Combined Slip Correction

`S_x` = Slip ratio

`alpha` = Slip angle

Normalize Slip factors

`S_x^star` = `(S_x)/(S__peak)`

`alpha^star` = `(alpha)/(alpha__peak)`

The maximum lateral force occurs at peak slip angle(`alpha`) or slip ratio (`S_x`)

The resultant slip is given my 

`S^star  = sqrt(((S_x^star)^2)+(alpha^2)))`

`S_x` modified = `S^star * S__peak`

`alpha` modified = `S^star * alpha__peak`

These equations were substituted in the above Pacjeka magic formula and determined

`F_x = F_xo  *  ((S_x^star)/(S^star))`

`F_y = F_yo  *  ((alpha^star)/(S^star))`

clear all
close all
clc

% Effect of skid on lateral force

Sxp = 0.09; % peak slip is 9%
alpha_p = 4; % peak slip angle is 4 degrees

% Pac formula coefficients
By = 0.27; Cy = 1.2; Dy = 2900; Ey = -1.6; Shy = 0; Svy = 0; %lateral force
Bx = 25; Cx = 1.15; Dx = 3200; Ex = -0.4; Shx = 0; Svx = 0; %longitudinal force

%iterate through four slip ratios (0,3,6,9)
for i = 1:4,% slip ratio loop
    Sx = (i-1)*0.03;
    for j = 1:21 % slip angle loop (o to 20 deg)
        slip_angle(j) = (j-1)+0.00001 ;
        alpha = slip_angle(j);
        
        %normalized slip factors - normalization done because the slip angle and slip ratio are in different units
        
        Sx_star = Sx/Sxp;
        alpha_star = alpha/alpha_p;
        
        %resultant slip
        S_star = (((Sx_star)^2 + (alpha_star)^2)^0.5);
        
        % modified slip factors
        Sx_mod = S_star*Sxp;
        alpha_mod = S_star*alpha_p;
        
        % Lateral force calculation using pac formula
        alpha_final = alpha_mod + Shy; % including horizantal shift
        fy(i,j) = (alpha_star/S_star)* Dy*sin(Cy*atan(By*alpha_final - Ey*(By*alpha_final - atan(By*alpha_final)))); % magic formula
        Fy(i,j) = fy(i,j) + Svy;
        
        
        % Longitudinal force calculation using pac formula
        Sx_final = Sx_mod + Shx; % including horizantal shift
        fx(i,j) = (Sx_star/S_star)*Dx*sin(Cx*atan(Bx*Sx_final - Ex*(Bx*Sx_final - atan(Bx*Sx_final))));% magic_formula
        Fx(i,j) = fx(i,j) + Svx; %including vertical shift
        
    end
end

Conclusions:

1. From the plots, it is clear that the lateral force reaches its peak value and later begins to saturate afterwards.

2. The maximum lateral force is achieved when the slip ratio is zero. And with an increase in slip ratio, the cornering capability reduces.

3. At zero skid ratio, the longitudinal force becomes zero and with an increase in skid ratio the longitudinal force increases.

 Tyre test data for curve fitting can be found here:

https://drive.google.com/file/d/0B2S2ZJKcURAmd0Vrdm1pSDhUY1k/view