MATLAB Examples

## Contents

```function OUTPUT = mfeval(PARAMETER_SOURCE,INPUTS,USE_MODE) %#codegen ```
```%MFEVAL evaluates Magic Formula 6.1.2 in steady state for series of input %variables. % % OUTPUT = mfeval(PARAMETER_SOURCE,INPUTS,USE_MODE) % % The formulation includes combined force/moment and turn slip calculation % ISO coordinate system is used in all calculations. % % PARAMETER_SOURCE refers to a MF-Tyre tyre property file (.TIR) containing % all the Magic Formula coefficients or to a structure with all the % parameters. % % INPUTS = [FZ KAPPA ALPHA GAMMA PHIT VX P* OMEGA*], where % FZ = normal load on the tyre [N] % KAPPA = longitudinal slip [dimensionless, -1: locked wheel] % ALPHA = side slip angle [rad] % GAMMA = inclination angle [rad] % PHIT = turn slip [1/m] % VX = forward velocity [m/s] % P* = pressure [Pa] % OMEGA* = rotational speed [rad/s] % % P* and OMEGA* are optional inputs. If they are not included pressure is % constant and equal to the nominal pressure on the TIR file and the % rotational speed is approximated. % % USE_MODE specifies the type of calculation performed: % 1: combined force/moment calculation % 2: combined force/moment calculation + turn slip % +10: revoke alpha_star definition % +20: include alpha_star definition % +100: include limit checks % +200: ignore limit checks % % For example: USE_MODE = 121 implies: % -combined force/moment % -include alpha_star % -include limit checks % % For normal situations turn slip may be neglected, meaning that the radius % of the path is close to infinity (straight line). % Alpha_star improves the accuracy of the model for very large slip angles % and possibly backward running of the wheel. % The limit checks verify that the inputs are inside the stable range of % the model and activates the low speed model. % % OUTPUT consists of 20 columns: % 1 - Fx: longitudinal force 11 - Vx: longitudinal velocity % 2 - Fy: lateral force 12 - P: pressure % 3 - Fz: normal force 13 - Re: effective rolling radius % 4 - Mx: overturning moment 14 - rho: tyre deflection % 5 - My: rolling resistance moment 15 - a: contact patch length % 6 - Mz: self aligning moment 16 - t: pneumatic trail % 7 - kappa: longitudinal slip 17 - longitudinal friction coefficient % 8 - alpha: side slip angle 18 - lateral friction coefficient % 9 - gamma: inclination angle 19 - omega: rot. speed % 10 - phit: turn slip 20 - Rl : loaded radius % % The equations here presented are published in the book: % Title: Tire and Vehicle Dynamics % Author: Hans Pacejka % Edition: 3, revised % Publisher: Elsevier, 2012 % ISBN: 0080970176, 9780080970172 % Length: 672 pages % % And in the following paper: % Besselink, I. J. M. , Schmeitz, A. J. C. and Pacejka, H. B.(2010) 'An % improved Magic Formula/Swift tyre model that can handle inflation % pressure changes', Vehicle System Dynamics, 48: 1, 337 — 352 % Link to this Article: DOI: 10.1080/00423111003748088 % URL: http://dx.doi.org/10.1080/00423111003748088 % % Script version: 1.4 % Date: 03/07/2017 % Code developed by: Marco Furlan % Email address 1: mfurlant@jaguarlandrover.com % Email address 2: marcofurlan92@gmail.com % % See the <a href="matlab:web(fullfile(mfevalRoot, 'html','index.html'))">documentation</a> for more details and examples ```

## Versions change log

```% Version | Date | Description % ------- ---- ----------- % V1.0 25/01/2017 -First implementation % % V1.1 31/01/2017 -Eqn (4.E40) Bt modified. % -Allow a structure of parameters to be % passed in instead of a TIR file name. % -Replace for loops in the coeff checks with % vector expressions. % -Rho calculation reviewed avoiding (NaNs) % in external effects. % -Coefficient checks can be disabled. % -Coefficients of friction are checked % before being exported. % % V1.2 16/02/2017 -Typo in the help text corrected. % -Mux and Muy warnings can be deactivated % with the coefficient checks. % -Rolling resistance (4.E70), pneumatic % trail (4.E73), Bt (4.E40), Dt (4.E43) and s % (4.E76) equations have been modified to % match TNO solver. % % V1.3 06/03/2017 -Code modified for code generation % -Simulink blocks with transients have been % included % -Use mode implementation modified % -Turn Slip equations modified % -Pneumatic trail equation fixed % % V1.4 03/07/2017 -First implementation of the low speed mode % -Publication of mfeval in MathWorks File % Exchange ```

## Future work

- Low speed model is not perfect, review the effect on pneumatic trail (t) and aligning torque (Mz) - Turn slip equations needs to be reviewed, specially for Mz

## Declare extrinsic functions

If mfeval is used in a Simulink block, the code generation software will generate code for the call to an extrinsic function. Functions 'warning' and 'loadTIR' are not supported for code generation. By adding coder.extrinsic('warning') we can to bypass code generation. Type : coder.screener('mfeval') into the command window to check

```coder.extrinsic('warning') coder.extrinsic('loadTIR') ```

## Check the arguments of the function

If the USE_MODE (3rd argument) is not specified, set it to 111

```if nargin < 3 USE_MODE = 111; warning('USE_MODE not defined, default mode 111 selected') end ```

## Apply Use Mode

Read the USE_MODE and apply the changes

```if isempty(USE_MODE) USE_MODE = 111; warning('USE_MODE is empty, default mode 111 selected') end % Check the number of digits of USE_MODE numdigits_USE_MODE = max(ceil(log10(abs(USE_MODE))),1); if numdigits_USE_MODE~=3 error('USE_MODE should contain exactly 3 digits') end % Save the elements of the string into separate variables hundreds = floor(USE_MODE/100); tens = floor((USE_MODE-hundreds*100)/10); units = floor((USE_MODE-hundreds*100-tens*10)); switch hundreds case 1 % 1: include limits checks useLimitsCheck = true; case 2 % 2: revoke limits checks useLimitsCheck = false; otherwise error('USE_MODE not contemplated') end switch tens case 1 % 1: revoke alpha_star definition useAlphaStar = false; case 2 % 2: include alpha_star definition useAlphaStar = true; otherwise error('USE_MODE not contemplated') end switch units case 1 % 1: combined force/moment calculation useTurnSlip = false; case 2 % 2: combined force/moment calculation + turn slip useTurnSlip = true; otherwise error('USE_MODE not contemplated') end ```

## Inputs (Parameters)

```if(ischar(PARAMETER_SOURCE)) % File location TIRparam = loadTIR(PARAMETER_SOURCE); % Call an external function to % read all the parameters of the TIR file. elseif(isstruct(PARAMETER_SOURCE)) TIRparam = PARAMETER_SOURCE; % Structure of parameter else error('PARAMETER_SOURCE not recognised') end % Low speed singularity epsilon = 1e-6; % [Eqn (4.E6a) Page 178 - Book] epsilonv = epsilon; epsilonx = epsilon; epsilonk = epsilon; epsilony = epsilon; epsilonr = epsilon; ```

## Inputs (Variables)

Unpack the input matrix into separate input variables

```Fz = INPUTS(:,1); % vertical load (N) kappa = INPUTS(:,2); % longitudinal slip (-) (-1 = locked wheel) alpha = INPUTS(:,3); % side slip angle (radians) gamma = INPUTS(:,4); % inclination angle (radians) phit = INPUTS(:,5); % turn slip (1/m) Vcx = INPUTS(:,6); % forward velocity (m/s) % Check the size of the inputs [mrows,ncolumns] = size(INPUTS); % If the pressure is specified by the user, grab it. % If is not, use the Inflation pressure of the TIR file if ncolumns > 6 p = INPUTS(:,7); % pressure (Pa) else p = ones(mrows,1).*TIRparam.INFLPRES; % pressure (Pa) end % IMPORTANT NOTE: Vx = Vcx [Eqn (7.4) Page 331 - Book] % It is assumed that the difference between the wheel centre longitudinal % velocity Vx and the longitudinal velocity Vcx of the contact centre is % negligible ```

## Limits

This code applies the appropriate limits to the input values to the model based on the MF6.1 limits specified in the TIR File

```% Save unlimited variables for other calculations Fz_unlimited = Fz; kappa_unlimited = kappa; if useLimitsCheck % Minimum Speed VXLOW = TIRparam.VXLOW; % Inflation Pressure Range PRESMIN = TIRparam.PRESMIN; PRESMAX = TIRparam.PRESMAX; % Vertical Force Range FZMIN = TIRparam.FZMIN; FZMAX = TIRparam.FZMAX; % Slip Angle Range ALPMIN = TIRparam.ALPMIN; ALPMAX = TIRparam.ALPMAX; % Inclination Angle Range CAMMIN = TIRparam.CAMMIN; CAMMAX = TIRparam.CAMMAX; % Long Slip Range KPUMIN = TIRparam.KPUMIN; KPUMAX = TIRparam.KPUMAX; % Speed limits (apply low slip mode) isLowSpeed = abs(Vcx) < VXLOW; % If Vcx is close to zero then SR is also 0 (linear digression) kappa(isLowSpeed) = kappa(isLowSpeed).*(Vcx(isLowSpeed)./VXLOW).*... sign(Vcx(isLowSpeed)); % Create a reduction factor for forces and moments to get zero output % at zero speed fraction = Vcx(isLowSpeed)./TIRparam.VXLOW; reduction_factor = sin(fraction.*1.5708); % 90 deg into rad isLowAndNegative = Vcx(isLowSpeed)<0; reduction_factor(isLowAndNegative) = -reduction_factor(isLowAndNegative); % Turn slip modifications phit = phit.*cos(alpha); % Empirically discovered % If the speed is negative, the turn slip is also negative isNegativeSpeed = Vcx<0; phit(isNegativeSpeed) = -phit(isNegativeSpeed); % If Vcx is close to zero, use linear digression phit(isLowSpeed) = phit(isLowSpeed).*abs(fraction); % Check Slip Angle limits isLowSlip = alpha < ALPMIN; alpha(isLowSlip) = ALPMIN; isHighSlip = alpha > ALPMAX; alpha(isHighSlip) = ALPMAX; % Check camber limits isLowCamber = gamma < CAMMIN; gamma(isLowCamber) = CAMMIN; isHighCamber = gamma > CAMMAX; gamma(isHighCamber) = CAMMAX; % Check Fz limits Fz(Fz < 0) = 0; isLowLoad = Fz < FZMIN & Fz > 0; Fz(isLowLoad) = FZMIN; isHighLoad = Fz > FZMAX & Fz > 0; Fz_unlimited = Fz; % Some formulas use the unsaturated Fz in TNO Fz(isHighLoad) = FZMAX; % Check pressure limits isLowPressure = p < PRESMIN; p(isLowPressure) = PRESMIN; isHighPressure = p > PRESMAX; p(isHighPressure) = PRESMAX; % Check slip ratio limits isLowSlipR = kappa < KPUMIN; kappa(isLowSlipR) = KPUMIN; isHighSlipR = kappa > KPUMAX; kappa(isHighSlipR) = KPUMAX; % Flag if anything is out of range. chk = isLowSlip | isHighSlip | isLowCamber | ... isHighCamber | isLowLoad | isHighLoad | isLowPressure | ... isHighPressure | isLowSlipR | isHighSlipR; if(any(chk)) warning (['Limits of the input values to the model exceeded.',... 'Values have been saturated.']); end if any(isLowSpeed) warning(['Speed input VX below the limit. ',... 'Low speed mode activated.']); end else % Not using limits checks % isLowSpeed = logical(zeros(length(Fz),1)); isLowSpeed = false(length(Fz),1); reduction_factor = 1; end ```

## First calculations

```% Basic calculations [Fzo_prime, dfz, dpi, alpha_star, gamma_star, alpha_prime, LMUX_star, LMUY_star, LMUX_prime, LMUY_prime] = calculateBasic(TIRparam,useAlphaStar, Fz, kappa, alpha, gamma, Vcx, p, epsilonv); % Effective rolling radius [Re, Romega, omega] = calculateRe(TIRparam, ncolumns, INPUTS, dpi, Fz_unlimited, kappa_unlimited); % Unlimited inputs used to match TNO results % Contact patch lengths [a, ~] = calculateContactPatch(TIRparam, Fz_unlimited, dpi); % Unlimited Fz used to match TNO results ```

## Turn slip calculations

Extension of the Model for Turn Slip, Page 183

```[zetaStr] = calculateZeta(TIRparam, useTurnSlip, useLimitsCheck, Fz, alpha, gamma, kappa, Vcx, omega, phit, alpha_star, gamma_star, Fzo_prime, dfz, dpi, LMUY_star, LMUY_prime, epsilonr, epsilonk); ```

## Fxo - Pure Longitudinal Force

(alpha = 0) Page 179 of the book

```[Fxo, mux, Kxk] = calculateFxo(TIRparam, useLimitsCheck, Fz, kappa, gamma, Vcx, dfz, dpi, epsilonx, LMUX_star, LMUX_prime, zetaStr, isLowSpeed, reduction_factor); ```

## Fyo - Pure Lateral Force

(kappa = 0) Page 179 of the book

```[Fyo, muy, Kya, Kygo, SHy, SVy, By, Cy] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, gamma_star, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr); ```

## Mzo - Pure Aligning Torque

(kappa = 0) Page 180 of the book

```[~, alphar, alphat, Dr, Cr, Br, Dt, Ct, Bt, Et, Kya_prime] = calculateMzo(TIRparam, useAlphaStar, useLimitsCheck, Fz, gamma, alpha_star, alpha_prime, gamma_star, Fzo_prime, dfz, dpi, Vcx, epsilonk, epsilony, Kya, Kygo, SHy, SVy, By, Cy, LMUY_star, LMUY_prime, zetaStr); ```

## Fx - Combined Longitudinal Force

(Combined Slip), Page 181 of the book

```[Fx] = calculateFx(TIRparam, useLimitsCheck, kappa, alpha_star, gamma_star, dfz, Fxo); ```

## Fy - Combined Lateral Force

(Combined Slip), Page 181 of the book

```[Fy, ~] = calculateFy(TIRparam, useLimitsCheck, Fz, alpha_star, kappa, gamma_star, dfz, Fyo, muy, zetaStr, isLowSpeed, reduction_factor); ```

## Mx - Overturning Couple

(Also see Section 4.3.5), Page 182 of the book and Page 12 of the paper

```[Mx] = calculateMx(TIRparam, Fz_unlimited, Fy, gamma, dpi); ```

## My - Rolling Resistance Moment

(See Eqns (9.236, 9.230, 9.231)), Page 182 of the book

```[My] = calculateMy(TIRparam, Fz_unlimited, Fx, Vcx, kappa, gamma, p, isLowSpeed, reduction_factor); ```

## Mz - Aligning Torque

(Combined slip), Page 182 of the book

```[Mz, t] = calculateMz(TIRparam, useAlphaStar, useLimitsCheck, alphar, alphat, Kxk, Kya_prime, kappa, Fz, Fy, Fx, Fzo_prime, dfz, gamma, Dr, Cr, Br, Dt, Ct, Bt, Et, alpha_star, Fyo, muy, alpha_prime, Vcx, dpi, epsilony, epsilonk, LMUY_star, LMUY_prime, isLowSpeed, reduction_factor); ```

Methodology explained by Henning Olsson (henning.olsson@calspan.com)

```[rho, Rl] = calculateRhoRl(TIRparam, omega, Fz_unlimited, Fx, Fy, dpi, Romega); % Unlimited Fz used to match TNO results ```

## Generate the output matrix

```% Check the sign of the coefficient of friction % The calculation of Fx and Fy is not affected by the sign of mux and muy if(useLimitsCheck) if any(mux < 0) mux = abs(mux); warning('Negative longitudinal coefficient of friction forced to be positive') end if any(muy < 0) muy = abs(muy); warning('Negative lateral coefficient of friction forced to be positive') end end % Preallocate OUTPUT variable OUTPUT = zeros(mrows,20); % Pack all the outputs OUTPUT(:,1) = Fx; OUTPUT(:,2) = Fy; OUTPUT(:,3) = Fz_unlimited; OUTPUT(:,4) = Mx; OUTPUT(:,5) = My; OUTPUT(:,6) = Mz; OUTPUT(:,7) = kappa; OUTPUT(:,8) = alpha; OUTPUT(:,9) = gamma; OUTPUT(:,10) = phit; OUTPUT(:,11) = Vcx; OUTPUT(:,12) = p; OUTPUT(:,13) = Re; OUTPUT(:,14) = rho; OUTPUT(:,15) = 2*a; OUTPUT(:,16) = t; OUTPUT(:,17) = mux; OUTPUT(:,18) = muy; OUTPUT(:,19) = omega; OUTPUT(:,20) = Rl; ```
```end ```

## Supplementary functions

Functions to calculate forces and moments

```function [Fzo_prime, dfz, dpi, alpha_star, gamma_star, alpha_prime, LMUX_star, LMUY_star, LMUX_prime, LMUY_prime] = calculateBasic(TIRparam,useAlphaStar, Fz, kappa, alpha, gamma, Vcx, p, epsilonv)%#codegen %[MODEL] LONGVL = TIRparam.LONGVL ; %Nominal speed %[OPERATING_CONDITIONS] NOMPRES = TIRparam.NOMPRES ; %Nominal tyre inflation pressure %[VERTICAL] FNOMIN = TIRparam.FNOMIN ; %Nominal wheel load %[SCALING_COEFFICIENTS] LFZO = TIRparam.LFZO ; % Scale factor of nominal (rated) load LMUX = TIRparam.LMUX ; % Scale factor of Fx peak friction coefficient LMUY = TIRparam.LMUY ; % Scale factor of Fy peak friction coefficient % New scaling factor in Pacejka 2012 with it's default value. This % scaling factor is not present in the standard MF6.1 TIR files. LMUV = 0 ; % Scale factor with slip speed Vs decaying friction % Rename the TIR file variables in the Pacejka style Fzo = FNOMIN; % Nominal (rated) wheel load pio = NOMPRES; % Reference pressure Vo = LONGVL; % Nominal speed (LONGVL) % Some basic calculations are done before starting to calculate forces and % moments. % (1) Velocities in point S (slip) Vsx = -kappa.*abs(Vcx); % [Eqn (4.E5) Page 181 - Book] Vsy = -tan(alpha).*abs(Vcx); % [Eqn (2.12) Page 67 - Book] and [(4.E3) Page 177 - Book] Vs = sqrt(Vsx.^2+Vsy.^2); % [Eqn (3.39) Page 102 - Book] -> Slip velocity of the slip point S % (2) Velocities in point C (contact) Vcy = Vsy; % Assumption from page 67 of the book, paragraph above Eqn (2.11) Vc = sqrt(Vcx.^2+Vcy.^2); % Velocity of the wheel contact centre C, Not described in the book but is the same as [Eqn (3.39) Page 102 - Book] % (3) Effect of having a tire with a different nominal load Fzo_prime = LFZO.*Fzo; % [Eqn (4.E1) Page 177 - Book] % (4) Normalized change in vertical load dfz = (Fz - Fzo_prime)./Fzo_prime; % [Eqn (4.E2a) Page 177 - Book] % (5) Normalized change in inflation pressure dpi = (p - pio)./pio; % [Eqn (4.E2b) Page 177 - Book] % (6) Use of alpha_star definition (Large slip angles correction) if useAlphaStar alpha_star = tan(alpha).*sign(Vcx); % [Eqn (4.E3) Page 177 - Book] else alpha_star = alpha; % To match TNO results end % (7) Spin due to the camber angle gamma_star = sin(gamma); % [Eqn (4.E4) Page 177 - Book] % (8) For the aligning torque at high slip angles Vc_prime = Vc + epsilonv.*sign(Vc); % [Eqn (4.E6a) Page 178 - Book] [sign(Vc) term explained on page 177] alpha_prime = acos(Vcx./Vc_prime); % [Eqn (4.E6) Page 177 - Book] % (9) Slippery surface with friction decaying with increasing (slip) speed LMUX_star = LMUX./(1 + LMUV.*Vs./Vo); % [Eqn (4.E7) Page 179 - Book] LMUY_star = LMUY./(1 + LMUV.*Vs./Vo); % [Eqn (4.E7) Page 179 - Book] % (10) Digressive friction factor % On Page 179 of the book is suggested Amu=10, but after comparing the use % of the scaling factors against TNO, Amu = 1 was giving perfect match Amu = 1; LMUX_prime = Amu.*LMUX_star./(1+(Amu-1).*LMUX_star); % [Eqn (4.E8) Page 179 - Book] LMUY_prime = Amu.*LMUY_star./(1+(Amu-1).*LMUY_star); % [Eqn (4.E8) Page 179 - Book] end function [zetaStr] = calculateZeta(TIRparam, useTurnSlip, useLimitsCheck, Fz, alpha, gamma, kappa, Vcx, omega, phit, alpha_star, gamma_star, Fzo_prime, dfz, dpi, LMUY_star, LMUY_prime, epsilonr, epsilonk)%#codegen % Declare extrinsic functions coder.extrinsic('warning') if useTurnSlip %[SCALING_COEFFICIENTS] LKY = TIRparam.LKY ; % Scale factor of Fy cornering stiffness LYKA = TIRparam.LYKA ; % Scale factor of alpha influence on Fx LKYC = TIRparam.LKYC ; % Scale factor of camber force stiffness LKZC = TIRparam.LKZC ; % %Scale factor of camber torque stiffness LMP = TIRparam.LMP ; % Scale factor of Parking Moment %[TURNSLIP_COEFFICIENTS] PDXP1 = TIRparam.PDXP1 ; %Peak Fx reduction due to spin parameter PDXP2 = TIRparam.PDXP2 ; %Peak Fx reduction due to spin with varying load parameter PDXP3 = TIRparam.PDXP3 ; %Peak Fx reduction due to spin with kappa parameter PKYP1 = TIRparam.PKYP1 ; %Cornering stiffness reduction due to spin PDYP1 = TIRparam.PDYP1 ; %Peak Fy reduction due to spin parameter PDYP2 = TIRparam.PDYP2 ; %Peak Fy reduction due to spin with varying load parameter PDYP3 = TIRparam.PDYP3 ; %Peak Fy reduction due to spin with alpha parameter PDYP4 = TIRparam.PDYP4 ; %Peak Fy reduction due to square root of spin parameter PHYP1 = TIRparam.PHYP1 ; %Fy-alpha curve lateral shift limitation PHYP2 = TIRparam.PHYP2 ; %Fy-alpha curve maximum lateral shift parameter PHYP3 = TIRparam.PHYP3 ; %Fy-alpha curve maximum lateral shift varying with load parameter PHYP4 = TIRparam.PHYP4 ; %Fy-alpha curve maximum lateral shift parameter PECP1 = TIRparam.PECP1 ; %Camber w.r.t. spin reduction factor parameter in camber stiffness PECP2 = TIRparam.PECP2 ; %Camber w.r.t. spin reduction factor varying with load parameter in camber stiffness QDTP1 = TIRparam.QDTP1 ; %Pneumatic trail reduction factor due to turn slip parameter QCRP1 = TIRparam.QCRP1 ; %Turning moment at constant turning and zero forward speed parameter QCRP2 = TIRparam.QCRP2 ; %Turn slip moment (at alpha=90deg) parameter for increase with spin QBRP1 = TIRparam.QBRP1 ; %Residual (spin) torque reduction factor parameter due to side slip QDRP1 = TIRparam.QDRP1 ; %Turn slip moment peak magnitude parameter if isfield(TIRparam,'QDRP2')==0 %Coefficient not found in the TIR file QDRP2 = 0 ; %Turn slip moment peak position parameter else QDRP2 = TIRparam.QDRP2 ; %Turn slip moment peak position parameter end %[LATERAL_COEFFICIENTS] PDY1 = TIRparam.PDY1 ; %Lateral friction Muy PDY2 = TIRparam.PDY2 ; %Variation of friction Muy with load PDY3 = TIRparam.PDY3 ; %Variation of friction Muy with squared camber PKY1 = TIRparam.PKY1 ; %Maximum value of stiffness Kfy./Fznom PKY2 = TIRparam.PKY2 ; %Load at which Kfy reaches maximum value PKY3 = TIRparam.PKY3 ; %Variation of Kfy./Fznom with camber PKY4 = TIRparam.PKY4 ; %Curvature of stiffness Kfy PKY5 = TIRparam.PKY5 ; %Peak stiffness variation with camber squared PKY6 = TIRparam.PKY6 ; %Fy camber stiffness factor PKY7 = TIRparam.PKY7 ; %Vertical load dependency of camber stiffness PVY3 = TIRparam.PVY3 ; %Variation of shift Svy./Fz with camber PVY4 = TIRparam.PVY4 ; %Variation of shift Svy./Fz with camber and load PPY1 = TIRparam.PPY1 ; %influence of inflation pressure on cornering stiffness PPY2 = TIRparam.PPY2 ; %influence of inflation pressure on dependency of nominal tyre load on cornering stiffness PPY3 = TIRparam.PPY3 ; %linear influence of inflation pressure on lateral peak friction PPY4 = TIRparam.PPY4 ; %quadratic influence of inflation pressure on lateral peak friction PPY5 = TIRparam.PPY5 ; %Influence of inflation pressure on camber stiffness RBY1 = TIRparam.RBY1 ; %Slope factor for combined Fy reduction RBY2 = TIRparam.RBY2 ; %Variation of slope Fy reduction with alpha RBY3 = TIRparam.RBY3 ; %Shift term for alpha in slope Fy reduction RBY4 = TIRparam.RBY4 ; %Influence of camber on stiffness of Fy combined RCY1 = TIRparam.RCY1 ; %Shape factor for combined Fy reduction REY1 = TIRparam.REY1 ; %Curvature factor of combined Fy REY2 = TIRparam.REY2 ; %Curvature factor of combined Fy with load RHY1 = TIRparam.RHY1 ; %Shift factor for combined Fy reduction RHY2 = TIRparam.RHY2 ; %Shift factor for combined Fy reduction with load %[ALIGNING_COEFFICIENTS] QDZ8 = TIRparam.QDZ8 ; %Variation of peak factor Dmr" with camber QDZ9 = TIRparam.QDZ9 ; %Variation of peak factor Dmr" with camber and load QDZ10 = TIRparam.QDZ10 ; %Variation of peak factor Dmr with camber squared QDZ11 = TIRparam.QDZ11 ; %Variation of Dmr with camber squared and load %[DIMENSION] Ro = TIRparam.UNLOADED_RADIUS ; %Free tyre radius % Speed limits (avoid zero speed) Vc_prime = Vcx; isLowSpeed = abs(Vcx) < TIRparam.VXLOW; signVcx = sign(Vcx(isLowSpeed)); signVcx(signVcx == 0) = 1; Vc_prime(isLowSpeed) = TIRparam.VXLOW * signVcx; % Singularity protected velocity, text on Page 184 psi_dot = phit.*Vc_prime; % Rearrange [Eqn (4.75) Page 183 - Book] % IMPORTANT NOTE: Eqn (4.75) has been modified % In chapter 2.2 "Definition of tire input quantities" in the Pacejka % book, it is assumed that the z-axis of the road coordinate system % "points downward normal to the road plane" (p. 62). Due to this % definition, Pacejka introduces the minus sign for the spin slip so % that a positive spin leads to a positive torque Mz (p. 68).But in % CarMaker (and other MBS software), the road coordinate system is % orientated differently. The z-axis points upward to the % road plane. Thus switching the signs is not needed here. epsilong = PECP1.*(1 + PECP2.*dfz); % [Eqn (4.90) Page 186 - Book] Camber reduction factor phi = (-1./Vc_prime).*(psi_dot-(1-epsilong).*omega.*sin(gamma)); % [Eqn (4.76) Page 184 - Book] zeta3 = cos(atan(PKYP1.*Ro.^2.*phi.^2)); % [Eqn (4.79) Page 185 - Book] Kygo = Fz.*(PKY6 + PKY7 .*dfz).*(1 + PPY5.*dpi).*LKYC; % (=dFyo./dgamma at alpha = gamma = 0) (= CFgamma) (4.E30) Kya = PKY1.*Fzo_prime.*(1 + PPY1.*dpi).*(1 - PKY3.*abs(gamma_star)).*sin(PKY4.*atan((Fz./Fzo_prime)./((PKY2+PKY5.*gamma_star.^2).*(1+PPY2.*dpi)))).*zeta3.*LKY; % (= ByCyDy = dFyo./dalphay at alphay = 0) (if gamma =0: =Kyao = CFa) (PKY4=2)(4.E15) Kyao = PKY1.*Fzo_prime.*(1 + PPY1.*dpi).*(1 - PKY3.*abs(0)).*sin(PKY4.*atan((Fz./Fzo_prime)./((PKY2+PKY5.*0.^2).*(1+PPY2.*dpi)))).*zeta3.*LKY; % IMPORTANT NOTE: Explanation of the above equation, Kyao % Kya0 is the cornering stiffness when the camber angle is zero % (gamma=0) which is again the product of the coefficients By, Cy and % Dy at zero camber angle. Information from Kaustub Ragunathan, email: % carmaker-service-uk@ipg-automotive.com signKya = sign(Kya); signKya(signKya == 0) = 1; % If [Kya = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN] in Eqn 4.E39 signKyao = sign(Kyao); signKyao(signKyao == 0) = 1; % If [Kyao = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN] Kya_prime = Kya + epsilonk.*signKya; % (4.E39) [sign(Kya) term explained on page 177] Kyao_prime = Kyao + epsilonk.*signKyao; % epsilonk is a small factor added to avoid the singularity condition during zero velocity (equation 308, CarMaker reference Manual). muy = (PDY1 + PDY2 .* dfz).*(1 + PPY3.*dpi + PPY4 .*dpi.^2).*(1 - PDY3.*gamma_star.^2).*LMUY_star; % (4.E23) Cyk = RCY1; % (4.E63) Byk = (RBY1 + RBY4.*gamma_star.^2).*cos(atan(RBY2.*(alpha_star - RBY3))).*LYKA; % (> 0) (4.E62) Eyk = REY1 + REY2.*dfz; % (<=1) (4.E64) SHyk = RHY1 + RHY2.*dfz; % (4.E65) Gyko = cos(Cyk.*atan(Byk.*SHyk - Eyk.*(Byk.*SHyk - atan(Byk.*SHyk)))); % (4.E60) Gyk_k = cos(Cyk.*atan(Byk.*kappa - Eyk.*(Byk.*kappa-atan(Byk.*kappa))))./Gyko; % (> 0)(4.E59) Byp = PDYP1.*(1 + PDYP2.*dfz).*cos(atan(PDYP3.*tan(alpha))); % [Eqn (4.78) Page 185 - Book] Bxp = PDXP1.*(1 + PDXP2.*dfz).*cos(atan(PDXP3.*kappa)); % [Eqn (4.106) Page 188 - Book] CHyp = PHYP1; % (>0) [Eqn (4.85) Page 186 - Book] DHyp = (PHYP2 + PHYP3.*dfz).*sign(Vcx); % [Eqn (4.86) Page 186 - Book] EHyp = PHYP4; % (<=1) [Eqn (4.87) Page 186 - Book] % Limits check if(useLimitsCheck) if EHyp > 1 EHyp=1; warning('Ex over limit (>1), Eqn(4.E14)') end end KyRpo = Kygo./(1-epsilong); % Eqn (4.89) BHyp = KyRpo./(CHyp.*DHyp.*Kyao_prime); %[Eqn (4.88) Page 186 - Book] SHyp = DHyp.*sin(CHyp .* atan(BHyp.*Ro.*phi - EHyp.*(BHyp.*Ro.*phi - atan(BHyp.*Ro.*phi)))).*sign(Vcx); % [Eqn (4.80) Page 185 - Book] Mzp_inf = QCRP1.*muy.*Ro.*Fz.*sqrt(Fz./Fzo_prime).*LMP; % [Eqn (4.95) Page 187 - Book] % Mzp_inf should be always > 0 negative_Mzp_inf = Mzp_inf < 0; Mzp_inf(negative_Mzp_inf) = 1e-6; Mzp90 = Mzp_inf.*(2./pi).*atan(QCRP2.*Ro.*abs(phit)).*Gyk_k; % [Eqn (4.103) Page 188 - Book] CDrp = QDRP1; % (>0) [Eqn (4.96) Page 187 - Book] EDrp = QDRP2; % (<=1) [Eqn (4.97) Page 187 - Book] DDrp = Mzp_inf./sin(0.5.*pi.*CDrp); % [Eqn (4.94) Page 187 - Book] Kzgro = Fz.*Ro.*(QDZ8.*QDZ9.*dfz + (QDZ10 + QDZ11.*dfz.*abs(gamma))).*LKZC; %[Eqn (4.99) Page 187 - Book] BDrp = Kzgro./(CDrp.*DDrp.*(1 - epsilong) + epsilonr.*sign(1 - epsilong)); % [Eqn (4.98) Page 187 - Book] [sign(1 - epsilong) term is explained on page 177] Drp = DDrp.*sin(CDrp.*atan(BDrp.*Ro.*(-phi) - EDrp.*(BDrp.*Ro.*(-phi) - atan(BDrp.*Ro.*(-phi))))); % Eqn (4.93) Page 187 - Book] % IMPORTANT NOTE: % Drp equation modified to operate with -phi, so that a positive spin % leads to a positive torque Mz (p. 68). Original equation below: % Drp = DDrp.*sin(CDrp.*atan(BDrp.*Ro.*(phi) - EDrp.*(BDrp.*Ro.*(phi) - atan(BDrp.*Ro.*(phi))))); % Eqn (4.93) Page 187 - Book] zeta0 = 0; % [Eqn (4.83) Page 186 - Book] zeta1 = cos(atan(Bxp.*Ro.*phi)); % [Eqn (4.105) Page 188 - Book] zeta2 = cos(atan(Byp.*(Ro.*abs(phi)+ PDYP4.*sqrt(Ro.*abs(phi))))); % [Eqn (4.77) Page 184 - Book] SVyg = Fz.*(PVY3 + PVY4.*dfz).*gamma_star.* LKYC .* LMUY_prime .* zeta2; % (4.E28) zeta4 = 1 + SHyp - SVyg./Kya_prime; % [Eqn (4.84) Page 186 - Book] zeta5 = cos(atan(QDTP1.*Ro.*phi)); %[Eqn (4.91) Page 186 - Book] zeta6 = cos(atan(QBRP1.*Ro.*phi)); % [Eqn (4.102) Page 188 - Book] signDrp = sign(Drp); signDrp(signDrp == 0) = 1; % If [Drp = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN] in Eqn 4.104 x = Mzp90./(-abs(Drp)+epsilonr.*signDrp); % Argument of Eqn (4.104) [sign(Drp) term is explained on page 177] % Text from Pacejka book page 188: (In any case, the argument of % (4.104) should remain <1.) This is because the inputs of "acos" % shpuld be between +1 to -1 to aboid complex solutions. if any(x > 1) || any(x < -1) isAbove = x>1; isBelow = x<-1; x(isAbove) = 1; x(isBelow) = -1; end zeta7 = (2./pi).*acos(x); % [Eqn (4.104) Page 188 - Book] % IMPORTANT NOTE: % zeta7 equation modified to operate with -(abs(Drp)) so that a % positive spin leads to a positive torque Mz (p. 68). Original % equation below: % zeta7 = (2./pi).*acos(Mzp90./(abs(Drp)+epsilonr.*signDrp)); % [Eqn (4.104) Page 188 - Book] zeta8 = 1 + Drp; else % No turn slip and small camber angles % First paragraph on page 178 of the book zeta0 = 1; zeta1 = 1; zeta2 = 1; zeta3 = 1; zeta4 = 1; zeta5 = 1; zeta6 = 1; zeta7 = 1; zeta8 = 1; end % Pack zeta in a structure zetaStr.zeta0 = zeta0; zetaStr.zeta1 = zeta1; zetaStr.zeta2 = zeta2; zetaStr.zeta3 = zeta3; zetaStr.zeta4 = zeta4; zetaStr.zeta5 = zeta5; zetaStr.zeta6 = zeta6; zetaStr.zeta7 = zeta7; zetaStr.zeta8 = zeta8; end function [Fxo, mux, Kxk] = calculateFxo(TIRparam, useLimitsCheck, Fz, kappa, gamma, Vcx, dfz, dpi, epsilonx, LMUX_star, LMUX_prime, zetaStr, isLowSpeed, reduction_factor)%#codegen % Declare extrinsic functions coder.extrinsic('warning') % Unpack zeta structure zeta1 = zetaStr.zeta1; %[SCALING_COEFFICIENTS] LCX = TIRparam.LCX ; % Scale factor of Fx shape factor LEX = TIRparam.LEX ; % Scale factor of Fx curvature factor LKX = TIRparam.LKX ; % Scale factor of Fx slip stiffness LHX = TIRparam.LHX ; % Scale factor of Fx horizontal shift LVX = TIRparam.LVX ; % Scale factor of Fx vertical shift %[LONGITUDINAL_COEFFICIENTS] PCX1 = TIRparam.PCX1 ; %Shape factor Cfx for longitudinal force PDX1 = TIRparam.PDX1 ; %Longitudinal friction Mux at Fznom PDX2 = TIRparam.PDX2 ; %Variation of friction Mux with load PDX3 = TIRparam.PDX3 ; %Variation of friction Mux with camber squared PEX1 = TIRparam.PEX1 ; %Longitudinal curvature Efx at Fznom PEX2 = TIRparam.PEX2 ; %Variation of curvature Efx with load PEX3 = TIRparam.PEX3 ; %Variation of curvature Efx with load squared PEX4 = TIRparam.PEX4 ; %Factor in curvature Efx while driving PKX1 = TIRparam.PKX1 ; %Longitudinal slip stiffness Kfx./Fz at Fznom PKX2 = TIRparam.PKX2 ; %Variation of slip stiffness Kfx./Fz with load PKX3 = TIRparam.PKX3 ; %Exponent in slip stiffness Kfx./Fz with load PHX1 = TIRparam.PHX1 ; %Horizontal shift Shx at Fznom PHX2 = TIRparam.PHX2 ; %Variation of shift Shx with load PVX1 = TIRparam.PVX1 ; %Vertical shift Svx./Fz at Fznom PVX2 = TIRparam.PVX2 ; %Variation of shift Svx./Fz with load PPX1 = TIRparam.PPX1 ; %linear influence of inflation pressure on longitudinal slip stiffness PPX2 = TIRparam.PPX2 ; %quadratic influence of inflation pressure on longitudinal slip stiffness PPX3 = TIRparam.PPX3 ; %linear influence of inflation pressure on peak longitudinal friction PPX4 = TIRparam.PPX4 ; %quadratic influence of inflation pressure on peak longitudinal friction Cx = PCX1.*LCX; % (> 0) (4.E11) mux = (PDX1 + PDX2.*dfz).*(1 + PPX3.*dpi + PPX4.*dpi.^2).*(1 - PDX3.*gamma.^2).*LMUX_star; % (4.E13) Dx = mux.*Fz.*zeta1; % (> 0) (4.E12) Kxk = Fz.*(PKX1 + PKX2.*dfz).*exp(PKX3.*dfz).*(1 + PPX1.*dpi + PPX2.*dpi.^2).*LKX; % (= BxCxDx = dFxo./dkx at kappax = 0) (= Cfk) (4.E15) signDx = sign(Dx); signDx(signDx == 0) = 1; % If [Dx = 0] then [sign(0) = 0]. This is done to avoid [Kxk / 0 = NaN] in Eqn 4.E16 Bx = Kxk./(Cx.*Dx + epsilonx.*signDx); % (4.E16) [sign(Dx) term explained on page 177] SHx = (PHX1 + PHX2.*dfz).*LHX; % (4.E17) SVx = Fz.*(PVX1 + PVX2.*dfz).*LVX.*LMUX_prime.*zeta1; % (4.E18) kappax = kappa + SHx; % (4.E10) Ex = (PEX1 + PEX2.*dfz + PEX3.*dfz.^2).*(1 - PEX4.*sign(kappax)).*LEX; % (<=1) (4.E14) % Limits check if(useLimitsCheck) if(any(Ex > 1)) warning('Ex over limit (>1), Eqn(4.E14)') Ex(Ex > 1) = 1; end end Fxo = Dx.*sin(Cx.*atan(Bx.*kappax-Ex.*(Bx.*kappax-atan(Bx.*kappax))))+SVx; % (4.E9) % Backward speed check Fxo(Vcx < 0) = -Fxo(Vcx < 0); % Low speed model if any(isLowSpeed) Fxo(isLowSpeed) = Fxo(isLowSpeed).* reduction_factor ; end end function [Fyo, muy, Kya, Kygo, SHy, SVy, By, Cy] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, gamma_star, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr)%#codegen % Declare extrinsic functions coder.extrinsic('warning') % Unpack zeta structure zeta0= zetaStr.zeta0; zeta2= zetaStr.zeta2; zeta3= zetaStr.zeta3; zeta4= zetaStr.zeta4; %[SCALING_COEFFICIENTS] LCY = TIRparam.LCY ; % Scale factor of Fy shape factor LEY = TIRparam.LEY ; % Scale factor of Fy curvature factor LKY = TIRparam.LKY ; % Scale factor of Fy cornering stiffness LHY = TIRparam.LHY ; % Scale factor of Fy horizontal shift LVY = TIRparam.LVY ; % Scale factor of Fy vertical shift LKYC = TIRparam.LKYC ; % Scale factor of camber force stiffness %[LATERAL_COEFFICIENTS] PCY1 = TIRparam.PCY1 ; %Shape factor Cfy for lateral forces PDY1 = TIRparam.PDY1 ; %Lateral friction Muy PDY2 = TIRparam.PDY2 ; %Variation of friction Muy with load PDY3 = TIRparam.PDY3 ; %Variation of friction Muy with squared camber PEY1 = TIRparam.PEY1 ; %Lateral curvature Efy at Fznom PEY2 = TIRparam.PEY2 ; %Variation of curvature Efy with load PEY3 = TIRparam.PEY3 ; %Zero order camber dependency of curvature Efy PEY4 = TIRparam.PEY4 ; %Variation of curvature Efy with camber PEY5 = TIRparam.PEY5 ; %Variation of curvature Efy with camber squared PKY1 = TIRparam.PKY1 ; %Maximum value of stiffness Kfy./Fznom PKY2 = TIRparam.PKY2 ; %Load at which Kfy reaches maximum value PKY3 = TIRparam.PKY3 ; %Variation of Kfy./Fznom with camber PKY4 = TIRparam.PKY4 ; %Curvature of stiffness Kfy PKY5 = TIRparam.PKY5 ; %Peak stiffness variation with camber squared PKY6 = TIRparam.PKY6 ; %Fy camber stiffness factor PKY7 = TIRparam.PKY7 ; %Vertical load dependency of camber stiffness PHY1 = TIRparam.PHY1 ; %Horizontal shift Shy at Fznom PHY2 = TIRparam.PHY2 ; %Variation of shift Shy with load PVY1 = TIRparam.PVY1 ; %Vertical shift in Svy./Fz at Fznom PVY2 = TIRparam.PVY2 ; %Variation of shift Svy./Fz with load PVY3 = TIRparam.PVY3 ; %Variation of shift Svy./Fz with camber PVY4 = TIRparam.PVY4 ; %Variation of shift Svy./Fz with camber and load PPY1 = TIRparam.PPY1 ; %influence of inflation pressure on cornering stiffness PPY2 = TIRparam.PPY2 ; %influence of inflation pressure on dependency of nominal tyre load on cornering stiffness PPY3 = TIRparam.PPY3 ; %linear influence of inflation pressure on lateral peak friction PPY4 = TIRparam.PPY4 ; %quadratic influence of inflation pressure on lateral peak friction PPY5 = TIRparam.PPY5 ; %Influence of inflation pressure on camber stiffness Kya = PKY1.*Fzo_prime.*(1 + PPY1.*dpi).*(1 - PKY3.*abs(gamma_star)).*sin(PKY4.*atan((Fz./Fzo_prime)./((PKY2+PKY5.*gamma_star.^2).*(1+PPY2.*dpi)))).*zeta3.*LKY; % (= ByCyDy = dFyo./dalphay at alphay = 0) (if gamma =0: =Kyao = CFa) (PKY4=2)(4.E25) SVyg = Fz.*(PVY3 + PVY4.*dfz).*gamma_star.* LKYC .* LMUY_prime .* zeta2; % (4.E28) Kygo = Fz.*(PKY6 + PKY7 .*dfz).*(1 + PPY5.*dpi).*LKYC; % (=dFyo./dgamma at alpha = gamma = 0) (= CFgamma) (4.E30) signKya = sign(Kya); signKya(signKya == 0) = 1; % If [Kya = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN] in Eqn 4.E27 SHy = (PHY1 + PHY2.*dfz).* LHY + ((Kygo .*gamma_star - SVyg)./(Kya + epsilonk.*signKya)).*zeta0 +zeta4 -1; % (4.E27) [sign(Kya) term explained on page 177] SVy = Fz.*(PVY1 + PVY2.*dfz).*LVY.*LMUY_prime.*zeta2 + SVyg; % (4.E29) alphay = alpha_star + SHy; % (4.E20) Cy = PCY1.*LCY; % (> 0) (4.E21) muy = (PDY1 + PDY2 .* dfz).*(1 + PPY3.*dpi + PPY4 .*dpi.^2).*(1 - PDY3.*gamma_star.^2).*LMUY_star; % (4.E23) Dy = muy.*Fz.*zeta2; % (4.E22) signAlphaY = sign(alphay); signAlphaY(signAlphaY == 0) = 1; Ey = (PEY1 + PEY2.*dfz).*(1 + PEY5.*gamma_star.^2 - (PEY3 + PEY4.*gamma_star).*signAlphaY).*LEY; % (<=1)(4.E24) % Limits check if(useLimitsCheck) if(any(Ey > 1)) warning('Ey over limit (>1), Eqn(4.E24)') Ey(Ey > 1) = 1; end end signDy = sign(Dy); signDy(signDy == 0) = 1; % If [Dy = 0] then [sign(0) = 0]. This is done to avoid [Kya / 0 = NaN] in Eqn 4.E26 By = Kya./(Cy.*Dy + epsilony.*signDy); % (4.E26) [sign(Dy) term explained on page 177] Fyo = Dy .* sin(Cy.*atan(By.*alphay-Ey.*(By.*alphay - atan(By.*alphay))))+ SVy; % (4.E19) % Backward speed check % Only for alpha_star if(useAlphaStar) Fyo(Vcx < 0) = -Fyo(Vcx < 0); end end function [Mzo, alphar, alphat, Dr, Cr, Br, Dt, Ct, Bt, Et, Kya_prime] = calculateMzo(TIRparam, useAlphaStar, useLimitsCheck, Fz, gamma, alpha_star, alpha_prime, gamma_star, Fzo_prime, dfz, dpi, Vcx, epsilonk, epsilony, Kya, Kygo, SHy, SVy, By, Cy, LMUY_star, LMUY_prime, zetaStr)%#codegen % Declare extrinsic functions coder.extrinsic('warning') % Unpack zeta structure zeta0= zetaStr.zeta0; zeta2= zetaStr.zeta2; zeta5= zetaStr.zeta5; zeta6= zetaStr.zeta6; zeta7= zetaStr.zeta7; zeta8= zetaStr.zeta8; %[DIMENSION] Ro = TIRparam.UNLOADED_RADIUS ; %Free tyre radius %[SCALING_COEFFICIENTS] LKY = TIRparam.LKY ; % Scale factor of Fy cornering stiffness LTR = TIRparam.LTR ; % Scale factor of peak of pneumatic trail LRES = TIRparam.LRES ; % Scale factor for offset of residual torque LKZC = TIRparam.LKZC ; % %Scale factor of camber torque stiffness %[ALIGNING_COEFFICIENTS] QBZ1 = TIRparam.QBZ1 ; %Trail slope factor for trail Bpt at Fznom QBZ2 = TIRparam.QBZ2 ; %Variation of slope Bpt with load QBZ3 = TIRparam.QBZ3 ; %Variation of slope Bpt with load squared QBZ4 = TIRparam.QBZ4 ; %Variation of slope Bpt with camber QBZ5 = TIRparam.QBZ5 ; %Variation of slope Bpt with absolute camber QBZ9 = TIRparam.QBZ9 ; %Factor for scaling factors of slope factor Br of Mzr QBZ10 = TIRparam.QBZ10 ; %Factor for dimensionless cornering stiffness of Br of Mzr QCZ1 = TIRparam.QCZ1 ; %Shape factor Cpt for pneumatic trail QDZ1 = TIRparam.QDZ1 ; %Peak trail Dpt = Dpt.*(Fz./Fznom.*R0) QDZ2 = TIRparam.QDZ2 ; %Variation of peak Dpt" with load QDZ3 = TIRparam.QDZ3 ; %Variation of peak Dpt" with camber QDZ4 = TIRparam.QDZ4 ; %Variation of peak Dpt" with camber squared QDZ6 = TIRparam.QDZ6 ; %Peak residual torque Dmr" = Dmr./(Fz.*R0) QDZ7 = TIRparam.QDZ7 ; %Variation of peak factor Dmr" with load QDZ8 = TIRparam.QDZ8 ; %Variation of peak factor Dmr" with camber QDZ9 = TIRparam.QDZ9 ; %Variation of peak factor Dmr" with camber and load QDZ10 = TIRparam.QDZ10 ; %Variation of peak factor Dmr with camber squared QDZ11 = TIRparam.QDZ11 ; %Variation of Dmr with camber squared and load QEZ1 = TIRparam.QEZ1 ; %Trail curvature Ept at Fznom QEZ2 = TIRparam.QEZ2 ; %Variation of curvature Ept with load QEZ3 = TIRparam.QEZ3 ; %Variation of curvature Ept with load squared QEZ4 = TIRparam.QEZ4 ; %Variation of curvature Ept with sign of Alpha-t QEZ5 = TIRparam.QEZ5 ; %Variation of Ept with camber and sign Alpha-t QHZ1 = TIRparam.QHZ1 ; %Trail horizontal shift Sht at Fznom QHZ2 = TIRparam.QHZ2 ; %Variation of shift Sht with load QHZ3 = TIRparam.QHZ3 ; %Variation of shift Sht with camber QHZ4 = TIRparam.QHZ4 ; %Variation of shift Sht with camber and load PPZ1 = TIRparam.PPZ1 ; %effect of inflation pressure on length of pneumatic trail PPZ2 = TIRparam.PPZ2 ; %Influence of inflation pressure on residual aligning torque SHt = QHZ1 + QHZ2.*dfz + (QHZ3 + QHZ4.*dfz).*gamma_star; % (4.E35) signKya = sign(Kya); signKya(signKya == 0) = 1; % If [Kya = 0] then [sign(0) = 0]. This is done to avoid [num / 0 = NaN] in Eqn 4.E38 Kya_prime = Kya + epsilonk.*signKya; % (4.E39) [sign(Kya) term explained on page 177] SHf = SHy + SVy./Kya_prime; % (4.E38) alphar = alpha_star + SHf; % = alphaf (4.E37) alphat = alpha_star + SHt; % (4.E34) Dto = Fz.*(Ro./Fzo_prime).*(QDZ1 + QDZ2.*dfz).*(1 - PPZ1.*dpi).* LTR.*sign(Vcx); % (4.E42) % Dt = Dto.*(1 + QDZ3.*abs(gamma_star) + QDZ4.*gamma_star.^2).*zeta5; % (4.E43) % % IMPORTANT NOTE: The above original equation (4.E43) was not matching the % TNO solver. The coefficient Dt affects the pneumatic trail (t) and the % self aligning torque (Mz). % It was observed that when negative inclination angles where used as % inputs, there was a discrepancy between the TNO solver and mfeval. % This difference comes from the term QDZ3, that in the original equation % is multiplied by abs(gamma_star). But in the paper the equation is % different and the abs() term is not written. Equation (A60) from the % paper resulted into a perfect match with TNO. % Keep in mind that the equations from the paper don't include turn slip % effects. The term zeta5 has been added although it doesn't appear in the % paper. % Paper definition: Dt = (QDZ1 + QDZ2.*dfz).*(1 - PPZ1.*dpi).* (1 + QDZ3.*gamma + QDZ4.*gamma.^2).*Fz.*(Ro./Fzo_prime).*LTR.*zeta5; % (A60) % Bt = (QBZ1 + QBZ2.*dfz + QBZ3.*dfz.^2).*(1 + QBZ5.*abs(gamma_star) + QBZ6.*gamma_star.^2).*LKY./LMUY_star; %(> 0)(4.E40) % % IMPORTANT NOTE: In the above original equation (4.E40) it is used the % parameter QBZ6, which doesn't exist in the standard TIR files. Also note % that on page 190 and 615 of the book a full set of parameters is given % and QBZ6 doesn't appear. % The equation has been replaced with equation (A58) from the paper. % Paper definition: Bt = (QBZ1 + QBZ2.*dfz + QBZ3.*dfz.^2).*(1 + QBZ4.*gamma + QBZ5.*abs(gamma) ).*LKY./LMUY_star; %(> 0) (A58) Ct = QCZ1; % (> 0) (4.E41) Et = (QEZ1 + QEZ2.*dfz + QEZ3.*dfz.^2).*(1+(QEZ4+QEZ5.*gamma_star).*(2./pi).*atan(Bt.*Ct.*alphat)); % (<=1) (4.E44) % Limits check if(useLimitsCheck) if(any(Et > 1)) warning('Et over limit (>1), Eqn(4.E44)') Et(Et > 1) = 1; end end to = Dt.*cos(Ct.*atan(Bt.*alphat - Et.*(Bt.*alphat - atan(Bt.*alphat)))).*cos(alpha_prime); %t(aplhat)(4.E33) % Create a structure of zeta with phit = 0 zetaStr_sub0.zeta0 = 1; zetaStr_sub0.zeta1 = 1; zetaStr_sub0.zeta2 = 1; zetaStr_sub0.zeta3 = 1; zetaStr_sub0.zeta4 = 1; zetaStr_sub0.zeta5 = 1; zetaStr_sub0.zeta6 = 1; zetaStr_sub0.zeta7 = 1; zetaStr_sub0.zeta8 = 1; % Evaluate Fyo with gamma = 0 and phit = 0 [Fyo_sub0,~,~,~,~,~,~,~] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, 0, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr_sub0); Mzo_prime = -to.*Fyo_sub0; % gamma=phi=0 (4.E32) Dr = Fz.*Ro.*((QDZ6 + QDZ7.*dfz).*LRES.*zeta2 + ((QDZ8 + QDZ9.*dfz).*(1 + PPZ2.*dpi) + (QDZ10 + QDZ11.*dfz).*abs(gamma_star)).*gamma_star.*LKZC.*zeta0).*LMUY_star.*sign(Vcx).*cos(alpha_star) + zeta8 - 1; % (4.E47) Br = (QBZ9.*LKY./LMUY_star + QBZ10.*By.*Cy).*zeta6; % preferred: qBz9 = 0 (4.E45) Cr = zeta7; % (4.E46) % Evaluate Kya with gamma = 0 [~,~,Kya_sub0,~,~,~,~,~] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, 0, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr); Kzao = Dto.*Kya_sub0; %#ok<NASGU> % Kya for gamma=0 (= -dMzo./dAlphay at alphay=gamma=0)(=CMa) (4.E48) Kzgo = Fz.*Ro.*(QDZ8 + QDZ9.*dfz).*(1 + PPZ2.*dpi).*LKZC.*LMUY_star - Dto.*Kygo; %#ok<NASGU> %(=dMzo./dgamma at alpha = gamma = 0) (= CMg)(4.E49) Mzro = Dr.*cos(Cr.*atan(Br.*alphar)).*cos(alpha_prime); % =Mzr(alphar)(4.E36) Mzo = Mzo_prime + Mzro; % (4.E31) end function [Fx] = calculateFx(TIRparam, useLimitsCheck, kappa, alpha_star, gamma_star, dfz, Fxo)%#codegen % Declare extrinsic functions coder.extrinsic('warning') %[SCALING_COEFFICIENTS] LXAL = TIRparam.LXAL ; % Scale factor of alpha influence on Fx %[LONGITUDINAL_COEFFICIENTS] RBX1 = TIRparam.RBX1 ; %Slope factor for combined slip Fx reduction RBX2 = TIRparam.RBX2 ; %Variation of slope Fx reduction with kappa RBX3 = TIRparam.RBX3 ; %Influence of camber on stiffness for Fx combined RCX1 = TIRparam.RCX1 ; %Shape factor for combined slip Fx reduction REX1 = TIRparam.REX1 ; %Curvature factor of combined Fx REX2 = TIRparam.REX2 ; %Curvature factor of combined Fx with load RHX1 = TIRparam.RHX1 ; %Shift factor for combined slip Fx reduction Cxa = RCX1; % (4.E55) Exa = REX1 + REX2.*dfz; % (<= 1) (4.E56) % Limits check if(useLimitsCheck) if(any(Exa > 1)) warning('Exa over limit (>1), Eqn(4.E56)') Exa(Exa > 1) = 1; end end SHxa = RHX1; % (4.E57) Bxa = (RBX1 + RBX3.*gamma_star.^2).*cos(atan(RBX2.*kappa)).*LXAL; % (> 0) (4.E54) alphas = alpha_star + SHxa; % (4.E53) Gxao = cos(Cxa.*atan(Bxa.*SHxa-Exa.*(Bxa.*SHxa-atan(Bxa.*SHxa)))); % (4.E52) Gxa = cos(Cxa.*atan(Bxa.*alphas-Exa.*(Bxa.*alphas-atan(Bxa.*alphas))))./Gxao; % (> 0)(4.E51) Fx = Gxa.*Fxo; % (4.E50) end function [Fy,Gyk] = calculateFy(TIRparam, useLimitsCheck, Fz, alpha_star, kappa, gamma_star, dfz, Fyo, muy, zetaStr, isLowSpeed, reduction_factor)%#codegen % Declare extrinsic functions coder.extrinsic('warning') % Unpack zeta structure zeta2= zetaStr.zeta2; %[SCALING_COEFFICIENTS] LYKA = TIRparam.LYKA ; % Scale factor of alpha influence on Fx LVYKA = TIRparam.LVYKA ; % Scale factor of kappa induced Fy %[LATERAL_COEFFICIENTS] RBY1 = TIRparam.RBY1 ; %Slope factor for combined Fy reduction RBY2 = TIRparam.RBY2 ; %Variation of slope Fy reduction with alpha RBY3 = TIRparam.RBY3 ; %Shift term for alpha in slope Fy reduction RBY4 = TIRparam.RBY4 ; %Influence of camber on stiffness of Fy combined RCY1 = TIRparam.RCY1 ; %Shape factor for combined Fy reduction REY1 = TIRparam.REY1 ; %Curvature factor of combined Fy REY2 = TIRparam.REY2 ; %Curvature factor of combined Fy with load RHY1 = TIRparam.RHY1 ; %Shift factor for combined Fy reduction RHY2 = TIRparam.RHY2 ; %Shift factor for combined Fy reduction with load RVY1 = TIRparam.RVY1 ; %Kappa induced side force Svyk./Muy.*Fz at Fznom RVY2 = TIRparam.RVY2 ; %Variation of Svyk./Muy.*Fz with load RVY3 = TIRparam.RVY3 ; %Variation of Svyk./Muy.*Fz with camber RVY4 = TIRparam.RVY4 ; %Variation of Svyk./Muy.*Fz with alpha RVY5 = TIRparam.RVY5 ; %Variation of Svyk./Muy.*Fz with kappa RVY6 = TIRparam.RVY6 ; %Variation of Svyk./Muy.*Fz with atan(kappa) DVyk = muy.*Fz.*(RVY1 + RVY2.*dfz + RVY3.*gamma_star).*cos(atan(RVY4.*alpha_star)).*zeta2; % (4.E67) SVyk = DVyk.*sin(RVY5.*atan(RVY6.*kappa)).*LVYKA; % (4.E66) SHyk = RHY1 + RHY2.*dfz; % (4.E65) Eyk = REY1 + REY2.*dfz; % (<=1) (4.E64) % Limits check if(useLimitsCheck) if(any(Eyk > 1)) warning('Eyk over limit (>1), Eqn(4.E64)') Eyk(Eyk > 1) = 1; end end Cyk = RCY1; % (4.E63) Byk = (RBY1 + RBY4.*gamma_star.^2).*cos(atan(RBY2.*(alpha_star - RBY3))).*LYKA; % (> 0) (4.E62) kappas = kappa + SHyk; % (4.E61) Gyko = cos(Cyk.*atan(Byk.*SHyk - Eyk.*(Byk.*SHyk - atan(Byk.*SHyk)))); % (4.E60) Gyk = cos(Cyk.*atan(Byk.*kappas - Eyk.*(Byk.*kappas-atan(Byk.*kappas))))./Gyko; % (> 0)(4.E59) Fy = Gyk.*Fyo + SVyk; % (4.E58) % Low speed model if any(isLowSpeed) Fy(isLowSpeed) = Fy(isLowSpeed).* reduction_factor ; end end function [Mx] = calculateMx(TIRparam, Fz, Fy, gamma, dpi)%#codegen %[VERTICAL] Fzo = TIRparam.FNOMIN ; %Nominal wheel load %[DIMENSION] Ro = TIRparam.UNLOADED_RADIUS ; %Free tyre radius %[SCALING_COEFFICIENTS] LVMX = TIRparam.LVMX ; % Scale factor of Mx vertical shift LMX = TIRparam.LMX ; % Scale factor of overturning couple %[OVERTURNING_COEFFICIENTS] QSX1 = TIRparam.QSX1 ; %Vertical shift of overturning moment QSX2 = TIRparam.QSX2 ; %Camber induced overturning couple QSX3 = TIRparam.QSX3 ; %Fy induced overturning couple QSX4 = TIRparam.QSX4 ; %Mixed load lateral force and camber on Mx QSX5 = TIRparam.QSX5 ; %Load effect on Mx with lateral force and camber QSX6 = TIRparam.QSX6 ; %B-factor of load with Mx QSX7 = TIRparam.QSX7 ; %Camber with load on Mx QSX8 = TIRparam.QSX8 ; %Lateral force with load on Mx QSX9 = TIRparam.QSX9 ; %B-factor of lateral force with load on Mx QSX10 = TIRparam.QSX10 ; %Vertical force with camber on Mx QSX11 = TIRparam.QSX11 ; %B-factor of vertical force with camber on Mx QSX12 = TIRparam.QSX12 ; %Camber squared induced overturning moment QSX13 = TIRparam.QSX13 ; %Lateral force induced overturning moment QSX14 = TIRparam.QSX14 ; %Lateral force induced overturning moment with camber PPMX1 = TIRparam.PPMX1 ; %Influence of inflation pressure on overturning moment % Mx = Ro.*Fz.*(QSX1.*LVMX - QSX2.*gamma.*(1 + PPMX1.*dpi) + QSX3.*((Fy)/Fzo)... % + QSX4.*cos(QSX5.*atan((QSX6.*(Fz./Fzo)).^2)).*sin(QSX7.*gamma + QSX8.*atan... % (QSX9.*((Fy)/Fzo))) + QSX10.*atan(QSX11.*(Fz./Fzo)).*gamma).*LMX; %(4.E69) % % IMPORTANT NOTE: The above original equation (4.E69) is not used % because is not matching the results of the official TNO solver. Also, in % the book definition and in the official paper, parameters QSX12 QSX13 and % QSX14 are not used % Instead of using equations (4.E69) from the book or (A47) from the % official paper, it has been used the equation (49) described in the draft % paper of Besselink (Not the official paper). This draft can be downloaded % from: % https://pure.tue.nl/ws/files/3139488/677330157969510.pdf % purl.tue.nl/677330157969510.pdf % Draft paper definition: Mx = Ro .* Fz .* LMX .* (QSX1 .* LVMX - QSX2 .* gamma .* (1 + PPMX1 .* dpi) - QSX12 .* gamma .* abs(gamma) + QSX3 .* Fy ./ Fzo + ... QSX4 .* cos(QSX5 .* atan((QSX6 .* Fz ./ Fzo) .^ 2)) .* sin(QSX7 .* gamma + QSX8 .* atan(QSX9 .* Fy ./ Fzo)) + ... QSX10 .* atan(QSX11 .* Fz ./ Fzo) .* gamma) + Ro .* Fy .* LMX .* (QSX13 + QSX14 .* abs(gamma)); % (49) end function [My] = calculateMy(TIRparam, Fz_unlimited, Fx, Vcx, kappa, gamma, p, isLowSpeed, reduction_factor)%#codegen %[OPERATING_CONDITIONS] pio = TIRparam.NOMPRES ; %Nominal tyre inflation pressure %[MODEL] Vo = TIRparam.LONGVL ; %Nominal speed %[VERTICAL] Fzo = TIRparam.FNOMIN ; %Nominal wheel load %[DIMENSION] Ro = TIRparam.UNLOADED_RADIUS ; %Free tyre radius %[SCALING_COEFFICIENTS] LMY = TIRparam.LMY ; % Scale factor of rolling resistance torque %[ROLLING_COEFFICIENTS] QSY1 = TIRparam.QSY1 ; %Rolling resistance torque coefficient QSY2 = TIRparam.QSY2 ; %Rolling resistance torque depending on Fx QSY3 = TIRparam.QSY3 ; %Rolling resistance torque depending on speed QSY4 = TIRparam.QSY4 ; %Rolling resistance torque depending on speed .^4 QSY5 = TIRparam.QSY5 ; %Rolling resistance torque depending on camber squared QSY6 = TIRparam.QSY6 ; %Rolling resistance torque depending on load and camber squared QSY7 = TIRparam.QSY7 ; %Rolling resistance torque coefficient load dependency QSY8 = TIRparam.QSY8 ; %Rolling resistance torque coefficient pressure dependency % My = Fz.Ro*(QSY1 + QSY2.*(Fx./Fzo) + QSY3.*abs(Vcx./Vo) + QSY4.*(Vcx./Vo).^4 ... % +(QSY5 + QSY6.*(Fz./Fzo)).*gamma.^2).*((Fz./Fzo).^QSY7.*(p./pio).^QSY8).*LMY.; %(4.E70) % % IMPORTANT NOTE: The above equation from the book (4.E70) is not used % because is not matching the results of the official TNO solver. % This equation gives a positive output of rolling resistance, and in the % ISO coordinate system, My should be negative. Furthermore, the equation % from the book has an error, multiplying all the equation by Fz instead of % Fzo (first term). % Because of the previous issues it has been used the equation (A48) from % the paper. % Paper definition: My = -Ro.*Fzo.*LMY.*(QSY1 + QSY2.*(Fx./Fzo) + QSY3.*abs(Vcx./Vo) + QSY4.*(Vcx./Vo).^4 ... +(QSY5 + QSY6.*(Fz_unlimited./Fzo)).*gamma.^2).*((Fz_unlimited./Fzo).^QSY7.*(p./pio).^QSY8); %(A48) % Backward speed check My(Vcx < 0) = -My(Vcx < 0); % Negative SR check My(kappa<-1) = -My(kappa<-1); % Low speed model if any(isLowSpeed) My(isLowSpeed) = My(isLowSpeed).*reduction_factor; end end function [Mz,t] = calculateMz(TIRparam, useAlphaStar, useLimitsCheck, alphar, alphat, Kxk, Kya_prime, kappa, Fz, Fy, Fx, Fzo_prime, dfz, gamma, Dr, Cr, Br, Dt, Ct, Bt, Et, alpha_star, Fyo, muy, alpha_prime, Vcx, dpi, epsilony, epsilonk, LMUY_star, LMUY_prime, isLowSpeed, reduction_factor)%#codegen %[VERTICAL] Fzo = TIRparam.FNOMIN ; %Nominal wheel load %[DIMENSION] Ro = TIRparam.UNLOADED_RADIUS ; %Free tyre radius %[SCALING_COEFFICIENTS] LS = TIRparam.LS ; % Scale factor of moment arm of Fx LFZO = TIRparam.LFZO ; % Scale factor of nominal (rated) load %[ALIGNING_COEFFICIENTS] SSZ1 = TIRparam.SSZ1 ; %Nominal value of s./R0: effect of Fx on Mz SSZ2 = TIRparam.SSZ2 ; %Variation of distance s./R0 with Fy./Fznom SSZ3 = TIRparam.SSZ3 ; %Variation of distance s./R0 with camber SSZ4 = TIRparam.SSZ4 ; %Variation of distance s./R0 with load and camber % alphar_eq = sqrt(alphar.^2+(Kxk./Kya_prime).^2.*kappa.^2).*sign(alphar); % (4.E78) % alphat_eq = sqrt(alphat.^2+(Kxk./Kya_prime).^2.*kappa.^2).*sign(alphat); % (4.E77) % s = Ro.*(SSZ1 + SSZ2.*(Fy./Fzo_prime) + (SSZ3 + SSZ4.*dfz).*gamma_star).*LS; % (4.E76) % IMPORTANT NOTE: The equations 4.E78 and 4.E77 are not used due to small % differences discovered at negative camber angles with the TNO solver. % Instead equations A54 and A55 from the paper are used. % % IMPORTANT NOTE: The coefficient "s" (Equation 4.E76) determines the % effect of Fx into Mz. The book uses "Fzo_prime" in the formulation, % but the paper uses "Fzo". The equation (A56) from the paper has a better % correlation with TNO. alphar_eq = atan(sqrt(tan(alphar).^2+(Kxk./Kya_prime).^2.*kappa.^2)).*sign(alphar); % (A54) alphat_eq = atan(sqrt(tan(alphat).^2+(Kxk./Kya_prime).^2.*kappa.^2)).*sign(alphat); % (A55) s = Ro.*(SSZ1 + SSZ2.*(Fy./Fzo) + (SSZ3 + SSZ4.*dfz).*gamma).*LS; % (A56) Mzr = Dr.*cos(Cr.*atan(Br.*alphar_eq)); % (4.E75) % Create a structure of zeta with phit = 0 zetaStr_sub0.zeta0 = 1; zetaStr_sub0.zeta1 = 1; zetaStr_sub0.zeta2 = 1; zetaStr_sub0.zeta3 = 1; zetaStr_sub0.zeta4 = 1; zetaStr_sub0.zeta5 = 1; zetaStr_sub0.zeta6 = 1; zetaStr_sub0.zeta7 = 1; zetaStr_sub0.zeta8 = 1; % Evaluate Gyk with gamma = 0 and phit = 0 [~,Gyk_sub0] = calculateFy(TIRparam,useLimitsCheck,Fz,alpha_star,kappa,0,dfz,Fyo,muy,zetaStr_sub0, isLowSpeed, reduction_factor); % Evaluate Fyo with gamma = 0 and phit = 0 [Fyo_sub0,~,~,~,~,~,~,~] = calculateFyo(TIRparam, useAlphaStar, useLimitsCheck, Fz, alpha_star, 0, Fzo_prime, dfz, dpi, Vcx, epsilony, epsilonk, LMUY_star, LMUY_prime, zetaStr_sub0); Fy_prime = Gyk_sub0.*Fyo_sub0; % (4.E74) t = Dt.*cos(Ct.*atan(Bt.*alphat_eq - Et.*(Bt.*alphat_eq - atan(Bt.*alphat_eq)))).*cos(alpha_prime); % (4.E73) t = t*LFZO; % IMPORTANT NOTE: the above equation is not written in any source, but "t" % is multiplied by LFZO in the TNO solver. This has been empirically % discovered. % Low speed model if any(isLowSpeed) t(isLowSpeed) = t(isLowSpeed).*reduction_factor; Mzr(isLowSpeed) = Mzr(isLowSpeed).*reduction_factor; end Mz_prime = -t.*Fy_prime; % (4.E72) Mz = Mz_prime + Mzr + s.*Fx; % (4.E71) end function [rho, Rl] = calculateRhoRl(TIRparam, omega, Fz_unlimited, Fx, Fy, dpi, Romega)%#codegen % Declare extrinsic functions coder.extrinsic('warning') %[MODEL] Vo = TIRparam.LONGVL ; %Nominal speed %[VERTICAL] Fzo = TIRparam.FNOMIN ; %Nominal wheel load Cz0 = TIRparam.VERTICAL_STIFFNESS ; %Tyre vertical stiffness Q_V2 = TIRparam.Q_V2 ; %Vertical stiffness increase with speed Q_FZ2 = TIRparam.Q_FZ2 ; %Quadratic term in load vs. deflection Q_FCX = TIRparam.Q_FCX ; %Longitudinal force influence on vertical stiffness Q_FCY = TIRparam.Q_FCY ; %Lateral force influence on vertical stiffness PFZ1 = TIRparam.PFZ1 ; %Pressure effect on vertical stiffness %[DIMENSION] Ro = TIRparam.UNLOADED_RADIUS ; %Free tyre radius % Model parameters as QFZ1 that normally aren't present in the TIR files Q_FZ1 = sqrt((Cz0.*Ro./Fzo).^2 - 4.* Q_FZ2); % Rearranging [Eqn (4) Page 2 - Paper] % Split Eqn (A3.3) Page 619 of the Book into different bits: speed_effect = Q_V2 .* (Ro ./ Vo) .*abs(omega); fx_effect = (Q_FCX .* Fx ./ Fzo).^2; fy_effect = (Q_FCY .* Fy ./ Fzo).^2; pressure_effect = (1+PFZ1.*dpi); % Joining all the effects except tyre deflection terms: external_effects = (1 + speed_effect - fx_effect - fy_effect) .* pressure_effect.*Fzo; % Equation (A3.3) can be written as: % Fz = (Q_FZ2*(rho/Ro)^2 + Q_FZ1*(rho/Ro)) * external_effects % Rearranging all the terms we end up with a quadratic equation as: % ax^2 + bx + c = 0 % Q_FZ2*(rho/Ro)^2 + Q_FZ1*(rho/Ro) -(Fz/(external_effects)) = 0 % Note: use of capital letters to avoid confusion with contact patch % lengths "a" and "b" A = Q_FZ2; B = Q_FZ1; C = -(Fz_unlimited ./ (external_effects)); if all((B^2 - 4*A*C) > 0) x = (-B + sqrt(B^2 - 4*A.*C)) ./ (2*A); else warning('No positive solution for rho calculation') x = (-B + sqrt(B^2 - 4*A.*C)) ./ (2*A); end rho = x.*Ro; % The loaded radius is the free-spinning radius minus the deflection Rl = Romega - rho; % Eqn A3.2 Page 619 - Book assuming positive rho at all the time end function [a, b] = calculateContactPatch(TIRparam, Fz_unlimited, dpi)%#codegen % Rename the TIR file variables in the Pacejka style Ro = TIRparam.UNLOADED_RADIUS; % Free tyre radius w = TIRparam.WIDTH; % Nominal width of the tyre Cz0 = TIRparam.VERTICAL_STIFFNESS; % Vertical stiffness PFZ1 = TIRparam.PFZ1 ; %Pressure effect on vertical stiffness % Nominal stiffness (pressure corrected) Cz = Cz0 .* (1 + PFZ1.*dpi); % [Eqn (5) Page 2 - Paper] - Vertical stiffness adapted for tyre inflation pressure %[CONTACT_PATCH] Q_RA1 = TIRparam.Q_RA1 ; %Square root term in contact length equation Q_RA2 = TIRparam.Q_RA2 ; %Linear term in contact length equation Q_RB1 = TIRparam.Q_RB1 ; %Root term in contact width equation Q_RB2 = TIRparam.Q_RB2 ; %Linear term in contact width equation % (12) Contact patch length a = Ro.*(Q_RA2.*(Fz_unlimited./(Cz.*Ro)) + Q_RA1.*sqrt(Fz_unlimited./(Cz.*Ro))); %[Eqn (9) Page 3 - Paper] b = w.*(Q_RB2.*(Fz_unlimited./(Cz.*Ro)) + Q_RB1.*(Fz_unlimited./(Cz.*Ro)).^(1/3)); %[Eqn (10) Page 3 - Paper] (Not used later in this code) end function [Re, Romega, omega] = calculateRe(TIRparam, ncolumns, INPUTS, dpi, Fz_unlimited, kappa_unlimited)%#codegen % Excerpt from OpenTIRE MF6.1 implementation % Date: 2016-12-01 % Prepared for Marco Furlan/JLR % Questions: henning.olsson@calspan.com % Forward velocity (m/s) Vcx = INPUTS(:,6); % Rename the TIR file variables in the Pacejka style Fzo = TIRparam.FNOMIN; % Nominal (rated) wheel load Ro = TIRparam.UNLOADED_RADIUS; % Free tyre radius Vo = TIRparam.LONGVL; % Nominal speed (LONGVL) Cz0 = TIRparam.VERTICAL_STIFFNESS; % Vertical stiffness PFZ1 = TIRparam.PFZ1 ; %Pressure effect on vertical stiffness BREFF = TIRparam.BREFF ; %Low load stiffness effective rolling radius DREFF = TIRparam.DREFF ; %Peak value of effective rolling radius FREFF = TIRparam.FREFF ; %High load stiffness effective rolling radius Q_RE0 = TIRparam.Q_RE0 ; %Ratio of free tyre radius with nominal tyre radius Q_V1 = TIRparam.Q_V1 ; %Tyre radius increase with speed % Nominal stiffness (pressure corrected) Cz = Cz0 .* (1 + PFZ1.*dpi); % [Eqn (5) Page 2 - Paper] - Vertical stiffness adapted for tyre inflation pressure % Check if omega is one of the inputs % If it is, use it to calculate Re, otherwise it can be approximated with a % short iteration if ncolumns > 7 % Omega is one of the inputs omega = INPUTS(:,8); % rotational speed (rad/s) Romega = Ro .* (Q_RE0 + Q_V1.*((omega.*Ro)./Vo).^2); % [Eqn (1) Page 2 - Paper] - Centrifugal growth of the free tyre radius % Eff. Roll. Radius Re = Romega -(Fzo./Cz) .* ( DREFF .*atan( BREFF.*(Fz_unlimited./Fzo)) + FREFF.*(Fz_unlimited./Fzo)); % [Eqn (7) Page 2 - Paper] else % Omega is not specified and is going to be approximated % Initial guess of Re based on something slightly less than R0 Re = (Ro .* 0.965); for i=1:10 % Use the most up to date Re to calculate an omega % omega = Vcx ./ Re; % Old definition of Henning without kappa, not valid for brake and drive omega = (kappa_unlimited.*Vcx+Vcx)./Re; % [Eqn (2.5) Page 65 - Book] % Then we calculate free-spinning radius Romega = Ro .* (Q_RE0 + Q_V1.*((omega.*Ro)./Vo).^2); % [Eqn (1) Page 2 - Paper] - Centrifugal growth of the free tyre radius % Effective Rolling Radius Re = Romega -(Fzo./Cz) .* ( DREFF .*atan( BREFF.*(Fz_unlimited./Fzo)) + FREFF.*(Fz_unlimited./Fzo)); % [Eqn (7) Page 2 - Paper] end end end ```