Pacjeka combined slip

                                    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)`

                                      Combined Slip Correction

`S_x` = Slip ratio

`alpha` = Slip angle

Normalize Slip factors

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

`alpha^*` = `(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^*  = sqrt(((S_x^*)^2)+(alpha^2)))`

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

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

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

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

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

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) = (i-1)*0.03;
    for j = 1:21 % slip angle loop (o to 20 deg)
        slip_angle(j) = (j-1)+0.0000001 ;
        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) = 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

clear all
close all
clc

Sxp = 0.09; % peak slip is 9%
alpha_p = 4; % peak slip angle is 4 degrees
By = 0.27; Cy = 1.27; Dy = 2900; Ey = -1.6; Shy = 0; Svy = 0;
Bx = 25; Cx = 1.15; Dx = 3200; Ex = -0.4; Shx = 0; Svx = 0;

% Iterate through four slip ratios

for i = 1:4 % slip ratio loop (0%, 3%, 6%, 9%)
    Sx = (i-1)*0.03;
    for j = 1:21 % slip angles loop 0 to 20 degrees
        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;
        
        % Calculating cornering force
        
        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; % including vertical shift
        
        % Calculating longitudinal force
        
        Sx_final = Sx_mod + Shy; % 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) + Svy; % including vertical shift
        
    end
end

figure(1);
plot(slip_angle, Fy, 'linewidth',2)
xlabel('Slip angle (deg)');
ylabel('Fy (N)');
legend('skid = 0','skid = 0.03', 'skid = 0.06', 'skid = 0.09');

figure(2);
plot(slip_angle, fx, 'linewidth',2);
xlabel('Slip angle (deg)');
ylabel('Fx (N)');
legend('skid = 0','skid = 0.03', 'skid = 0.06', 'skid = 0.09');

% Part 2  Effect of Slip angle on longitudinal force

clear all
Sxp = 0.09; % peak slip is 9%
alpha_p = 4; % peak slip angle is 4 degrees
By = 0.27; Cy = 1.27; Dy = 2900; Ey = -1.6; Shy = 0; Svy = 0;
Bx = 25; Cx = 1.15; Dx = 3200; Ex = -0.4; Shx = 0; Svx = 0;

for i = 1:4
    slip_angle(i) = (i-1)*2; % slip angles (0 deg, 2 deg, 4 deg, 6 deg);
    alpha = slip_angle(i);
    for j = 1:30
        Sx(j) = (j-1)*0.02+0.00001;              
        Sx_star = Sx(j)/Sxp; % normalized slip factors
        alpha_star = alpha/alpha_p; % normalized slip factors
        S_star = ((Sx_star)^2+(alpha_star)^2)^0.5; % Resultant slip
        Sx_mod = S_star*Sxp; % modified slip factors
        alpha_mod = S_star*alpha_p; % modified slip factors
        
        % calculating cornering forces
        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; % including vertical shift
        
        % calculating longitudinal force
        
        Sx_final = Sx_mod + Shy; % 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)+ Svy; % including vertical shift
    end
end

figure(3);
plot(Sx*100, Fx,'linewidth',2);
xlabel('Skid (%)');
ylabel('Fx (N)')
legend('alpha = 0','alpha = 2 deg','alpha = 4 deg','alpha = 6 deg');

figure(4);
plot(Sx*100, Fy,'linewidth',2);
xlabel('Skid (%)');
ylabel('Fy (N)')
legend('alpha = 0','alpha = 2 deg','alpha = 4 deg','alpha = 6 deg');