Lack of precision trying to calculate IRI using a quarter car model

Hello there, I'm trying to calculate IRI (international roughness index) using a quarter car model. As you can see at the end of the code has 2 differents methods to get IRI, the first one is "IRI" wich uses the vertical displacement and the second one uses vertical velocity. To compare I'm using Proval Software and sometimes the second method gives me a number close to proval, but sometimes not and i don't know why. Any help will be welcome guys!!

*Note: the "IRI3" variable i'm trying to filter the output signal just to see what's happen

clear all
close all

%Restrições que não podem ser alteradas no código
%V = 22,22 m/s
%kt = rigidez pneu, ms = massa suspensa, mus = massa ñ suspensa, 
%cs =coeficiente de amortecimento, ks = rigidez da mola

%k/ms = 63.3, kt/ms = 653, c/ms = 6, mus/ms = 0,15
%Parametros encontrados em P 53 little book of profiling
  
% Caracteristicas massa suspensa

k = 6330;                    % Rigidez da mola (N/m) entre massa suspensa e não suspensa
c = 600;                     % Fator de amortecimento massa suspensa  
ms = 100;                    % 1/4 da massa suspensa  (kg)

% Caracteristicas massa não suspensa

mus = 15;                    % 1/4 da massa não suspensa  (kg)
kus = 65300;                 % Rigidez do pneu  (N/m)



v = 80/3.6;                   % Velocidade do veiculo  (m/s)
  


                        
%Representação estado de espaço verificar "modern control engineering ogata
%p88

Aqcar = [0 1 0 0;-k/ms -c/ms k/ms c/ms;0 0 0 1;k/mus c/mus -(k+kus)/mus -c/mus];
Bqcar = [0 0 0 kus/mus]'; Cqcar = [1 0 0 0; 0 0 1 0; 0 1 0 0;0 0 0 1]; 
Dqcar = 0;
qcar = ss(Aqcar,Bqcar,Cqcar,Dqcar);

%Cqcar=[x1,v1,x2,v2] baseado na formulação acima

% Respota y(:,1) deslocamento da massa suspensa
% Respota y(:,2) deslocamento da massa não suspensa
% Respota y(:,3) deslocamento relativo entre massa suspensa e não suspensa
% (adicionar um absoluto para eliminar valores negativos)
% Respota y(:,4) velocidade relativa entre massa suspensa e não suspensa

% FORMULAÇÃO ESTÁ CORRETA ATÉ AQUI

x0 = [0 0 0 0];                         % estado inicial

load sg_esq_inteira                   % Perfil de pista determinado

dx = dx_danilo(2) - dx_danilo(1);             % Diferencial de espaço dx entre coordenadas dx=0.01

dt2 = dx/v;                             % Diferencial de tempo dt entre coordenadas dt2=0.01/22.22
                                        

tempo = 0:dt2:(length(dx_danilo)-1)*dt2;



% Respostas do modelo proposto, baseados na entrada Cqcar
 y = lsim(qcar,Sf_elevacao_esq,tempo,x0);                                                                       

 

 
% PLOT RESPOSTA

    figure(1);clf                            
    plot(tempo,y(:,2), 'b');  hold on 
    plot(tempo,y(:,1),'g')
    legend('deslocamento massa n suspensa','deslocamento massa suspensa')
    xlabel('tempo (s)')
    ylabel('deslocamento (m)')
    
    %%% AREA DE CALCULO DO IRI 

% suspenso=y(:,1)            nsuspenso=y(:,2) *L40,41


delta=sum(abs(y(:,1)-y(:,2))); %Soma dos deslocamentos absolutos descrito pela subtração dos deslocamentos da massa suspensa e da massa não suspensa


IRI=(1/(v*tempo(end)*1e-3))*delta*dx %%UNIDADE m/km[

delta2=sum(abs(y(:,3)-y(:,4))); %Soma das velocidades absolutas


IRI2=(1/(v*tempo(end)*1e-3))*delta2*dt2 %%UNIDADE m/km

a=1;
b=[1/4 1/4 1/4 1/4];
c=(y(:,3));
d=(y(:,4));

filtro= filter(b,a,c);
filtro2= filter(b,a,d);

delta3=sum(abs(filtro-filtro2)); %Soma das velocidades absolutas


IRI3=(1/(v*tempo(end)*1e-3))*delta3*dt2 %%UNIDADE m/km

[Surendra Chowdari Gatiganti]

Hello Lucas,

Can you share the logitudianl profile file that you are using. so that i can run using my code that I have to better answer your question.

Thank you, have a good day.

-Suri

[Lucas Kffuri da Rocha]

Hello Surendra,

Sure I can share the profile. My code returns IRI by displacement method equals 10.2827 km/m and with the velocity method 5.4306

Proval software gives me a 3.664 km/m and Sayers method 3.65 km/m

dx_danilo is the dx variable and Sf_elevacao_esq the elevation variable

[Surendra Chowdari Gatiganti]

Hello Lucas,

Sounds like a silly question. But, is your Sf_elevacao_esq smoothend with moving average method with 250mm base width? I can see in your code that you are not smoothening it, but the Proval software does before simulating quater car.

Thank you.

[Lucas Kffuri da Rocha]

Hello Surendra, isn't a silly quest at all. To be honest, I'm confused about when to use 250 mm filter. In that case, for example, doesn't make any difference in proval either turn on or turn off the filter, the result is the same 3.664 km/m. If I'm not mistaken, in "guidelines for conducting and calibrating" has a recommendation to use filter when the dx is shorter than 0.25m

[Surendra Chowdari Gatiganti]

I think you might want to read this paper "On the Calculation of International Roughness Index from Longitudinal Road Profile"

[Surendra Chowdari Gatiganti]

[e-mail] send me an email on this, if you want to discuss about this. Thanks

Answer 1

[Mehmet Rizelioglu]

Hi Lucas, could you find answer your question? i have a problem like yours. can you help me to solve my problem?

[Lucas Kffuri da Rocha]

Hello Mehmet, yes I found the answer. You can find the final code in my paper: https://repositorio.utfpr.edu.br/jspui/handle/1/28063

If remains any doubts about it, you can send me an email [e-mail]