Implementation of the 1987 Magic Formula Tyre Model in Matlab

Magic Formula Tyre Model

The Magic Formula Tyre Model is an empirical tyre model, comprising of mathematical equations that can closely fit tyre test rig measurements. Despite, the tyre model does not take theoretical considerations into the tyre’s structures nor materials, it calculates the tyre forces and moments with reasonable accuracy at minimal calculation cost, therefore, making it appealing for use in simulation studies.

In this page, one of the early version of the Magic Formula Tyre Model is introduced and implemented in Matlab, considering pure slip cases.

1987 Magic Formula

The Magic Formula tyre model made its first appearance in 1987 by Bakker et al. [1]. The model is applicable to the longitudinal and lateral tyre forces \(F_x, \; F_y\) and the aligning moment \( M_z \). It has the capacity to consider combined slip conditions where both the longitudinal slip \( \kappa\) and the slip angle \(\alpha\) exist at the same time.

For the pure slip conditions, where only the longitudinal or lateral force exist,

Lateral force \( F_y\) \[ F_y = D \sin(C \arctan(B\phi) ) + \Delta S_v \] where \[ \left\{ \begin{array}{l} \begin{align} \phi &= (1-E)(\alpha + \Delta S_h ) + \dfrac{E}{B}\arctan(B (\alpha + \Delta S_h) )\\ D &= a_1 F_z^2 + a_2 F_z \\ C &= 1.30 \\ BCD &= a_3 \sin (a_4 \arctan(a_5 F_z)) (1 – a_{12} |\gamma|) \\ E &= a_6 F_z^2 + a_7 F_z + a_8 \\ \Delta S_h &= a_9 \gamma \\ \Delta S_v &= (a_{10} F_z^2 + a_{11} F_z) \gamma \end{align} \end{array} \right. \]

Brake force F_x \[ F_x = D \sin(C \arctan(B\phi) ) \] where \[ \left\{ \begin{array}{l} \begin{align} \phi &= (1-E)\kappa + \dfrac{E}{B}\arctan(B\kappa) \\ D &= b_1 F_z^2 + b_2 F_z \\ C &= 1.65 \\ BCD &= \dfrac{b_3 F_z^2 + b_4 F_z}{\mathrm{e}^{b_5 F_z}} \\ E &= b_6 F_z^2 + b_7 F_z + b_8 \\ \end{align} \end{array} \right. \]

The \( a \) and \(b\) parameters for the tyre studied by Bakker et al. [1] are quoted below.

Definitions of slip quantities in use

The Adapted SAE tyre axis system [2] is used, where the direction of the positive slip angle \(\alpha\) is the opposite of the SAE system. Consequently, the cornering stiffness \(C > 0 \) and \( \alpha > 0 \) produces positive lateral force \(F_y >0 \).

A tyre is in free rolling when it is pulled with a small force just enough to overcome the rolling resistance. For a free rolling tyre, the effective rolling radius \(r_e\) is defined as \[ r_e = \dfrac{V_x}{\Omega_o} \] where \(V_x\) is the longitudinal component of the velocity vector and \(\Omega_o\) corresponding angular velocity in free rolling.
When a torque is applied to the tyre, the longitudinal slip is defined as \[ \kappa = \; – \dfrac{V_x – \Omega r_e}{V_x} = \; – \dfrac{\Omega_o – \Omega}{\Omega_o} \] \(\kappa > 0\) corresponds with positive longitudinal force \(F_x > 0\), i.e. tractive force.
\( \kappa < 0 \) corresponds with negative longitudinal force \(F_x < 0 \), i.e. brake force.
If the wheel lock up \( \kappa = -1 \) since \( \Omega = 0 \).

Matlab function implementation

For the lateral force \(F_y\), the following Matlab function was written.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

function [Fy] = MF87Fy(alpha, Fz, gamma, Fyparameters)
% Fz - vertical load in [N]
% gamma - camber angle in [deg]
% alpha - slip angle in [deg]
 
% Fyparameters - .mat file that contains parameters e.g. 'MF87FyParameters.mat'
 
%initialise parameters a1 to a11
load(Fyparameters);
 
%unit conversion of inputs
Fz_kN = Fz./1000;
 
%calculation of coefficients
Dy = a1 * Fz_kN^2 + a2 * Fz_kN;
Cy = 1.30;
BCDy = a3 * sin(a4 *atan(a5*Fz_kN) ).*(1-a12*abs(gamma));
By = BCDy / (Cy*Dy);
Ey = a6 * Fz_kN^2 + a7 * Fz_kN + a8;
Shy = a9*gamma;
Svy = (a10*Fz_kN^2 + a11*Fz_kN)*gamma;
 
ByPhiy = By*(1-Ey)*(alpha+Shy) + Ey*atan(By*(alpha+Shy));
 
Fy = Dy*sin(Cy*atan(ByPhiy)) + Svy;
 
end

Similarly, for the longitudinal force \( F_x\),

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

function [Fx] = MF87Fx(kappa, Fz, Fxparameters)
% Fz - vertical load in [N]
% kappa - longitudinal slip in [%]
 
% Fxparameters - .mat file that contains parameters e.g. 'MF87FxParameters.mat'
 
%initialise parameters a1 to a11
load(Fxparameters);
 
%unit conversion of inputs
Fz_kN = Fz./1000;
 
%calculation of coefficients
Dx = b1 * Fz_kN^2 + b2 * Fz_kN;
Cx = 1.65;
BCDx = (b3 * Fz_kN^2 + b4 * Fz_kN)*exp(b5*Fz_kN*(-1));
Bx = BCDx / (Cx*Dx);
Ex = b6 * Fz_kN^2 + b7 * Fz_kN + b8;
 
BxPhix = Bx*(1-Ex)*kappa + Ex*atan(Bx*kappa);
 
Fx = Dx*sin(Cx*atan(BxPhix));
 
end

Results

By calling the functions for multiple times at sweeping longitudinal slip \( \kappa \) or slip angle \(\alpha \), the pure slip tyre force curves can be drawn.
The following figures are the plotted tyre force curves.

References

[1] Bakker et al. Tyre Modelling for Use in Vehicle Dynamics Studies, 1987.
[2] H. Pacejka, Tire and Vehicle Dynamics, 2012.