İç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 10100 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önder