Random road profile into ode45

11 次查看(过去 30 天)
Donghun Lee
Donghun Lee 2020-8-3
编辑: Drishti about 18 hours 前
clc,clear all
k_s = 26400;
m = 483; %Mass
f_n = sqrt(k_s/m)/(2*pi); %Natural frequency in Hz
%% Road profile
% spatial frequency (n0) cycles per meter
Omega0 = 0.1; %%%%conventional value of spatial frequency(n0)?
% psd ISO (used for formula 8)
Gd_0 = 32 * (10^-6);
% waveviness
w = 2;
% road length
L = 250;
%delta n
N = 1000;
Omega_L = 0.004;
Omega_U = 4;
delta_n = 1/L; % delta_n = (Omega_U - Omega_L)/(N-1);
% spatial frequency band
Omega = Omega_L:delta_n:Omega_U;
%PSD of road
Gd = Gd_0.*(Omega./Omega0).^(-w);
% calculate amplitude using formula(8) in the article
%Amp = sqrt(2*Gd*delta_n); %%%from Eq. 7?
%calculate amplitude using simplified formula(9) in the article
k = 3; %%%upper limit A and lower limit B k=3?
%Amp = sqrt(delta_n) * (2^k) * (10^-3) * (Omega0./Omega);
Amp = sqrt(delta_n) * (2^k) * (10^-3) * (Omega0./Omega);
%random phases
Psi = 2*pi*rand(size(Omega));
% x abicsa from 0 to L
x1 = 0:0.25:250;
% road sinal
h= zeros(size(x1));
plot(x1, h*1000);
%% ode45
T = 110;
x0 = [0,0];
for iv=1:length(x1)
h(iv) = sum( Amp.*cos(2*pi*Omega*x1(iv) + Psi) );
f = @(t,x) [ x(2); -( k_s*(x(1)-1000*h(iv))/ m ) ];
[t, x] = ode45(f,[100,T],x0);
end
plot(t,x(:,1));
Hi, I generated a random road profile (h) and it gives such a figure as below,
As you can see, there are quite lots of elevation in this road profile. However, if I apply this into my ode45, it gives a graph as below,
this seems to be a sinosoidal graph, but I think this graph should not be like this because the road has lots of elevation during the journey.
Can any one tell me what is wrong with my code please?

回答(1 个)

Drishti
Drishti about 18 hours 前
编辑:Drishti about 18 hours 前
Hi Donghun,
I understand that you are trying to plot a road profile followed by analyzing the system's response.
While reproducing the provided code on my end in MATLAB R2024b, I found that the potential reason for the system's response as a sinusoidal wave is because the 'ode45' is called inside the 'for' loop, recalculating the system's response at each point in 'x1'.
As a result, the simulation restarts at every step, whereas it should be integrated over the entire time span.
For better understanding, you can refer to the code snippet below:
for iv = 1:length(x1)
h(iv) = sum(Amp .* cos(2*pi*Omega*x1(iv) + Psi));
end
plot(x1, h*1000);
xlabel('Distance (m)');
ylabel('Elevation (mm)');
title('Road Profile');
%% ODE45
T = 110;
x0 = [0, 0];
f = @(t, x) [x(2); -(k_s * (x(1) - 1000 * interp1(x1, h, t)) / m)];
[t, x] = ode45(f, [0, T], x0);
For reference purpose, I have attached the output graph of system's response below :
I hope this resolves the issue.

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by