11 Haziran 2018 Pazartesi

Bisiklet Hareket Analizi – Uygulama (Matlab)


İçindekiler:

Arazi Eğiminin Belirlenmesi

Arazi yükseltisi konuma bağlı olarak biliniyor olsun. Arazi üzerindeki konum ile iz düşüm konumun yaklaşık olarak eşit olduğu belirtilmişti.
Konuma bağlı olarak bilinen yükselti değerleri bir dosyadan okunacak şekilde ayarlayalım ve dosyanın içeriği şu şekilde olsun:
0   10
15  11
30  12
45  11
60  9
Verilen sayılardan ilki s değeri, ikincisi ise h değeridir. Yukarıdaki değerleri yükseltiler.txt şeklinde çalışma klasörümüze kaydedelim.
%% dosya ayarları
fileName = 'yükseltiler.txt';  %yükseltilerin yer aldığı txt dosyası
a = [ 2 Inf ]; %bkz: fscanf()
spec = '%f\t%f'; %verilerin şablonu
%% dosya okuma
F = fopen(fileName, 'r');
DATA = fscanf(F, spec, a)';
fclose(F);
h_veri = DATA(:,2)'; %yükselti değerlerim
s_veri = DATA(:,1)'; %konum değerlerim
clear DATA F fileName A m n;
%% eğri uydurma
s = linspace(s_veri(1), s_veri(end));
H_s = createFit(s_veri, h_veri); %interpolasyonla tamamlanan H_s
h = H_s(s); %H_s şeklinde uydurulan fonksiyonda ...
            % ... s değerlerine karşılık h elde ediyorum

%% eğim
DH_s = sayisalTurev(s, h); %dh/ds in sayısal değerleri
theta = atan(DH_s); %theta'nın hesaplanması
%% grafik
subplot(2,1,1);
plot(s_veri, h_veri, 'o', s,h);
    title('Arazi');
    xlabel('Mesafe [m]');
    ylabel('Yükselti [m]');
subplot(2,1,2);
plot(s, theta*180/pi);
    title('Eğim');
    xlabel('Mesafe [m]');
    ylabel('Eğim [derece]');

Script dosyasında kullandığım iki fonksiyonun tanımlanması gerekiyor. createFit()ve sayisalTurev()fonksiyonlarının tanımları şu şekilde:
function [fitresult] = createFit(X, Y)
[xData, yData] = prepareCurveData( X, Y );
% Set up fittype and options.
ft = 'splineinterp';
% Fit model to data.
[fitresult] = fit( xData, yData, ft, 'Normalize', 'on' );

CreateFit fonksiyonu CurveFitting uygulamasıyla oluşturdum.

function [ Dy ] = sayisalTurev(X, Y)
n_x = length(X);
n_y = length(Y);
alta = 3;

if( n_x == n_y )
    Dy = zeros([1 n_x]);
    h = X(2) - X(1);
    for j = 1:n_x-alta;
       Dy(j) = ( Y(j+3)/3 - Y(j+2)*3/2 + Y(j+1)*3 - Y(j)*11/6 )/h;
    end
   
    for j = n_x-alta+1:n_x;
       Dy(j) = ( -Y(j-3)/3 + Y(j-2)*3/2 - Y(j-1)*3 + Y(j)*11/6 )/h;
    end
else
    Dy = -1;
    fprintf('Örten fonksiyon değil. \nX:%d\tY:%d\n', n_x, n_y);
end

end

Sayısal türev ileri ve geri fark polinomlarının türevini kullanan bir algoritma.  Türev verilerinin de uzunluk olarak örtüşmesi için hem ileri fark hem de geri fark kullanıldı. Nitekim tek birinin kullanılmasıyla veri adedi düşüyor. Scriptin çalıştırılmasıyla elde edilen grafik:

 

Arazide yol alan bisikletin hızı

 ve  biliniyor ve  bulunmak isteniyor. Bunun için birinci mertebeden diferansiyel denklem çözülmesi gerektiği belirtilmişti. Konuma bağlı hız değerleri incelendiğinden çözümün var olabilmesi için bisikletin geri gitmemesi, harmonik harekete başlamaması gerekir. Yani bisikletin aşamayacağı bir rampa sonucunda, durduğu noktadaki s değerinden itibaren bir çözüm yoktur.
Diferansiyel denklemini Düzeltilmiş Euler metoduyla çözen bir scripte yer vereceğim. Ardından aynı değerler için Matlab’ın sunduğu fonksiyonlardan kullanıp tekrar çıktı alacağım.
0   10
100    10
200    10
300    9
400    9

Yükselti değerlerimi yukarıdaki gibi değiştiriyorum.
%% dosya ayarları
fileName = 'yükseltiler.txt';  %yükseltilerin yer aldığı txt dosyası
a = [ 2 Inf ]; %bkz: fscanf()
spec = '%f\t%f'; %verilerin şablonu
%% dosya okuma
F = fopen(fileName, 'r');
DATA = fscanf(F, spec, a)';
fclose(F);
h_veri = DATA(:,2)'; %yükselti değerlerim
s_veri = DATA(:,1)'; %konum değerlerim
clear DATA F fileName A m n;
%% eğri uydurma
s = linspace(s_veri(1), s_veri(end));
H_s = createFit(s_veri, h_veri); %interpolasyonla tamamlanan H_s
h = H_s(s); %H_s şeklinde uydurulan fonksiyonda ...
            % ... s değerlerine karşılık h elde ediyorum
%% eğim
DH_s = sayisalTurev(s, h); %dh/ds in sayısal değerleri
theta = atan(DH_s); %theta'nın hesaplanması
%% değerler
g = 9.81; %yer çekimi ivmesi
b = 0.05; %sürtünme kuvveti katsayısı
M = 75; % toplam kütle
M2 = 70; % tekerlek kütlesi eksiltilmiş kütle
V_s = 1; % başlangıç hızı (m/s)
r_T = 0.36; %tekerlek yarıçapı
%% hız
theta_s = createFit(s, theta); %theta için fonksiyon oluşturalım
T = @(s) 1.7*(1+exp(-0.5*s)); %tekerleğe uygulanan tork
f = @(s,v) (T(s)/r_T-M*g*sin(theta_s(s)')-b*v-0.5)/M2/v; %dv/ds=f(s,v) denklemim
[sV, V] = dedd1(f,0.1,(s(end)-s(1))/0.1,0,V_s); %diferansiyel denklemin çözülmesi
%% grafik
subplot(3,1,1);
plot(s_veri, h_veri, 'o', s,h);
    title('Arazi');
    xlabel('Mesafe [m]');
    ylabel('Yükselti [m]');
subplot(3,1,2);
plot(sV, V*3.6);
    title('Sürat');
    xlabel('Mesafe [m]');
    ylabel('Sürat [km/h]');
subplot(3,1,3);
plot(s, T(s));
    title('Tekerleğe Verilen Tork');
    xlabel('Mesafe [m]');
    ylabel('Tork [Nm]');

Değerler kısmına kadar aynı dosyamı kullanabilirim. Tekerleğe uygulanan torkun konuma bağlı değişimi fonksiyon yardımı ile tanımlandı. Sabit bir tork değeri rampa için yetersiz, negatif eğim içinse fazla olduğundan dolayı başlangıç anında yüksek ardından sönüm yapmış bir fonksiyon olarak tanımlandı.  dedd1() fonksiyonu, düzeltilmiş euler yöntemi ile denklemi çözen fonksiyon.
function [ x, y ] = dedd1(fun,h,n,x0,y0)
% dy/dx = f(x,y) şeklinde olan fonksiyonlar için
% fun = @(x,y) f(x,y) şeklinde tanımlanmalıdır.
x = zeros([1 n]); x(1) = x0;
y = zeros([1 n]); y(1) = y0;
for i = 1:n-1;
    x1 = x(i)+h;
    d0 = fun(x(i), y(i));
    y1p = y(i)+h*d0;
    y1c = y(i)+h/2*(d0+fun(x1, y1p));
  
    x(i+1) = x1;
    y(i+1) = y1c;
end
clear x1 y1p d0 y1c i;
end

Scriptin çıktısı aşağıdaki gibidir.
Sadece diferansiyel denklemin çözümünde Matlab fonksiyonları kullanılırsa,
%% hız
theta_s = createFit(s, theta); %theta için fonksiyon oluşturalım
T = @(s) 4*(1+2*exp(-0.05*s))+100*theta_s(s)'; %tekerleğe uygulanan tork
f = @(s,v) (T(s)/r_T-M*g*sin(theta_s(s))-b*v-0.5)/M2/v; %dv/ds=f(s,v) denklemim
[sV, V] = ode23(f, [s(1) s(end)], V_s); %diferansiyel denklemin çözülmesi

Aynı çözüm (daha kısa sürede J) elde edilebilir.
Daha çok eğimli araziler için tork fonksiyonu düzeltilebilir.

 

Sürtünme kuvveti katsayısının belirlenmesi

Sürtünme kuvveti katsayısının belirlenmesi ile ilgili kısım teorik yazıda açıklanmıştı.
Bunun için düz yolda GPS kaydı yapılması gerekiyor. Bu kaydın aşağıdaki gibi olduğu veya aşağıdaki gibi düzenlendiği varsayılıyor.
40.20187586 29.08233377 172.8965066 0.42410147
40.2019745  29.08220868 171.2383162 0.2702889
40.20199417 29.08216902 168.4961052 0.47309476
40.20203649 29.08205445 164.0722264 0.48611397

İlk iki sütun konumun, 3. sütun yüksekliğin, 4. sütun ise hızın bilgisini içeriyor. Her satır 1’er saniyelik aralarla kaydedilmiş. Elbette yükseklikteki değişime bakılırsa düz yol değil J.
%% dosya
fileName = 'sürtünmeHesabı.txt'; %verilerin olduğu dosya
a = [ 4 Inf ];
hIndis = 3;
vIndis = 4;
spec = '%f\t%f\t%f\t%f'; %şablon
%% incelenecek aralık
timeConcept = [250 300]; %incelenecek olan satırlar
parts = 1000; %örnekleme sayısı
%% dosya oku
F = fopen(fileName, 'r');
DATA = fscanf(F, spec, a)';
m = size(DATA);
fclose(F);
h = DATA(:,hIndis)';
v = DATA(:,vIndis)';
if timeConcept(2)==0
   timeConcept = [1 m];
end
clear DATA F fileName A m n;
%% veriler
thisTime = timeConcept(1):timeConcept(2);
v_0 = v(thisTime); %incelenecek olan kısmın alınması
h_0 = h(thisTime);
a_0 = sayisalTurev(thisTime, v_0); %ivme
Dh_0 = sayisalTurev(thisTime, h_0); %dh/dt
s_0 = sayisalIntegral(thisTime, v_0); % s = s(t)
%% eğri uydurma
S_t = createFit(thisTime, s_0); % s = s(t) fonsiyonu interpolasyonu
H_t = createFit(thisTime, h_0); % h = h(t)
DH_t = createFit(thisTime, Dh_0); % h'(t)
V_t = createFit(thisTime, v_0); % v = v(t)
A_t = createFit(thisTime, a_0); % a = v'(t)
%% nümerik değerler
T = linspace(timeConcept(1), timeConcept(2), parts);
V = V_t(T); %incelenecek aralıktaki nümerik değerler
S = S_t(T);
DH = DH_t(T);
H = H_t(T);
A = A_t(T);
%% eğim
theta = atan(DH./V);
%% grafik
close all;
subplot(2,2,1);plot(S, H);
    title('Arazi');xlabel('Mesafe [m]');ylabel('Yükselti [m]');
subplot(2,2,3);plot(S, theta*180/pi);
    title('Eğim');xlabel('Mesafe [m]'); ylabel('Eğim [derece]');
subplot(2,2,2);plot(S, V);
    title('Hız-Mesafe');xlabel('Mesafe [m]'); ylabel('Sürat [m/s]');
subplot(2,2,4);plot(S, A);
    title('İvme-Mesafe');xlabel('Mesafe [m]'); ylabel('İvme [m/s^2]');
f = figure;
plot(S, A./V*70+9.81*sin(theta)./V);
xlabel('Mesafe [m]');ylabel('B')

Programla eğimin olmadığı ve tahrikin olmadığı aralık seçilerek b değeri elde edilebilir. Zaman değişkeninden Konum değişkenine geçiş zamana bağlı olan hızın integrali ile gerçekleştirilmiştir. Sayısal integral fonksiyonu aşağıdaki gibidir:

function [ F ] = sayisalIntegral(X, Y)
n_x = length(X);
n_y = length(Y);
if( n_x == n_y )
    F = zeros([1 n_x]);
    h = X(2) - X(1);
    F(1) = 0;
    for i = 1:n_x-1;
        F(i+1) = F(i)+Y(i)*h;
    end
else
    F = -1;
    fprintf('Örten fonksiyon değil. \nX:%d\tY:%d\n', n_x, n_y);
end
end

Program çıktısı:
Elle tutulur bir sabit b seçilen aralıkta yoktur. Daha hassas bir çalışma ve daha genel bir sürtünme modeli ile ilgili katsayılar script ile elde edilebilir.

Hiç yorum yok:

Yorum Gönderme