Odwrócone wahadło.

Układ składa się z wózka i zamocowanego na elastycznym przegubie pionowego wahadła.

DANE M - masa wózka 0.5 kg m - masa wahadła 0.5 kg b - współczynnik tarcia 0.1 N/m/sec l - połowa długoPci wahadla0.3 m I - moment bezwładnoPci wahadła 0.006 kg*m^2 F - siła przyłożona do wózka x - położenie wózka - kąt odchylenia wahadła z pozycji pionowej

Rozklad sił:

Równania ruchu:

1. dla wózka w kierunku poziomym

2.dla wahadła w kierunku poziomym

3. suma sił w kierunku prostopadłym do wahadła        

4. suma momentów względem środka wahadła

z 1 i 2 otrzymujemy:

z 2 i 3

Przykład 1

Napisać funkcję obliczającą kwadrat z dowolnej liczby.

 

f_kw.m

function a=f_kw(x)

a=x.^2

 

>> y=f_kw(10)

y=100

 

Przykład 2

Napisać funkcję obliczającą średnią arytmetyczną i geometryczną ciągu zapisanego w wektorze: a=[a1 a2 ... an]

 

srednia.m

function [sa,sg]=srednia(a)

n=length(a);

sa=sum(a)/n;

sg=(prod(a))^(1/n);

 

Możemy sprawdzić działanie tej funkcji

>> [a,g]=średnia([1:2:11])

gdzie wynikiem będzie liczba a (średnia arytmetyczna) i g (średnia geometryczna).

1. 1.        RÓWNANIE RÓŻNICZKOWE

 

Większość obiektów dynamicznych da się opisać równaniem różniczkowym postaci:

 

a_ny⁽ⁿ⁾+a_n-1y^(n-1)+... +a₁y¹+a₀y=u(t)

 

 

 

 

Każde równanie różniczkowe n-tego rzędu da się sprowadzić do układu n równań rzędu 1-go. Zróbmy następujące podstawienia:

y(t)=x₁(t)

x'₁=x₂

x'₂=x₃

:

x'_n-1=x_n

wówczas:

Otrzymamy układ równań 1-go rzędu następującej postaci:

który w notacji macierzowjma następującą postać:

=

 

Takie równanie macierzowe zapiszemy teraz w postaci funkcji Matlaba:

 

row_rozn.m

 

function dx=row_rozn(t,x)

global A B

dx=A*x+wymuszenie(t);

 

Funkcja reprezentująca wymuszenie (wektor sygnałów wejściowych będzie miała następującą postać:

 

wymuszenie.m

 

function u=wymuszenie(t)

if t >= 0

u = 1

else

u = 0

end

 

Aby numerycznie rozwiązać powyższe równanie różniczkowe w Matlabie (tzn. otrzymać wartości funkcji y dla odpowiadających im wartości chwil czasowych t należy użyć funkcji ode23 lub ode45. Skrypt rozwiązujący takie równanie powinien mieć następującą postać:

global A B

A=[ ]

B=[ ]

[t,y] = ode23(`rów_różn',[t₀ t_k],y₀)

gdzie: t,y - rozwiązania, wektory odpowiednio chwil czasowych i wartości funkcji,

`rów_różn' - nazwa funkcji definiującej rozwiązywane równanie różniczkowe,

t₀, t_k - przedział czasu, w którym ma być rozwiązywane równanie różniczkowe,

y₀ - warunki początkowe (wartość zmiennych x₁ …x_n w chwili t₀.

 

Przykład:

Oscylator liniowy z tłumieniem (obciążnik na sprężynie)

(przyjęto, iż nie występuje wymuszenie u=0)

 

Funkcja reprezentująca powyższy układ równań różniczkowych:

oscylator.m

 

function DX=oscylator(t,X)

global A

DX=A*X;

 

Skrypt rozwiązujący takie równanie:

 

global A ;

k=1;

m=1;

ၡ=0.5;

A=[0 1 ; -k/m -ၡ/m];

y0=[5 0]; (warunek początkowy dla położenia różny od 0, x₁=5)

[t,y]=ode23('oscylator',[0 20], y0);

plot(t,y(:,1));

Przykład:

Weźmy pod uwagę oscylator liniowy z tłumieniem (obciążnik na sprężynie) opisany następującym równaniem różniczkowym:

Przykład:

Jeżeli
, to

» l=[2 4] <Enter>

» m = [3 0 4 5] <Enter>

 

Ekran MATLABA powinien wyglądać tak:

» l=[2 4]

l =

2 4

» m = [3 0 4 5]

m =

3 0 4 5

 

W przestrzeni roboczej MATLABA istnieją teraz zmienne l (licznik transmitancji) i m (mianownik).

Wprowadź polecenie printsys(l,m) , otrzymasz

num/den =

2 s + 4

------------------

3 s^3 + 4 s + 5

 

W Matlabie istnieje również możliwość generacji standardowych modeli obiektów dynamicznych:

• -         obiekt oscylacyjny - ord2

• -         opóźnienie transportowe: pade(T₀,n) ,

gdzie: n - jest rzędem wielomianu, którym jest to opóźnienie przybliżone.

[l,m]=pade(T₀,n)

[A,B,C,D]=pade(T₀,n)

• -         model generowany losowo : rmodel(n) ; gdzie: n - jest to rząd transmitancji.

Charakterystyki czasowe:

•         Skokowa - polecenie step

step(l,m) lub step(A,B,C,D)

 

•         Impulsowa - polecenie impulse

impulse(l,m) lub impulse(A,B,C,D)

•       Dowolna - polecenie lsim

y=lsim(l,m,u,t)

y=lsim(A,B,C,D,u,t,x0)

gdzie: u(t) - wektor wartości sygnału wejściowego dla kolejnych chwil t

x0 - wektor wartości początkowe

 

Charakterystyki częstotliwościowe:

•         Bode'go - polecenie bode

bode(l,m)

•         Nyquist'a - polecenie nyquist

nyquist(l,m)

Przykład 1

Wyznaczyć odpowiedź skokową układu o transmitancji

dla ၷ_n = 1 oraz ၸ = 0.5. Wyniki symulacji wyświetlić na oscyloskopie.

Zadanie to zrealizowane zostanie w następującym układzie.

Rys. Kompletny model badanego układu

Przykład 2

Wyznaczyć odpowiedź skokową układu o transmitancji

dla ၷ_n = 2 oraz ၸ = 0.25. Czas trwania symulacji 10 sekund. Wyniki symulacji wyświetlić na oscyloskopie oraz zapisać do przestrzeni roboczej MATLABA i następnie:

 uzyskane wyniki zapisać do pliku dyskowego pod nazwą odp_skokowa

 zamknąć okno SIMULINKA

 wyczyścić przestrzeń roboczą MATLABA poleceniem clear

 odczytać umieszczone w pliku dyskowym dane i przedstawić je na wykresie

 uzyskany wykres zamieścić w dokumencie Worda

Zadanie to zrealizowane zostanie w następującym układzie.

Rys. Model układu z przykładu 2

Przykład 3

Wyznaczyć odpowiedź skokową układu o transmitancji

dla ၷ_n = 2 oraz ၸ = = 0.25, 0.5, 0.75. Czas trwania symulacji 10 sekund. Wyniki symulacji przedstawić na jednym wykresie.

Do realizacji tak postawionego zadania wykorzystany zostanie schemat pokazany na poprzednim rysunku i przechowywany pod nazwą uklad_IIrz.mdl. Skrypt wykonujący to zadanie jest następujący:

close all % Zamknięcie wszystkich okien graficznych

clear % Wyczyszczenie pamieci roboczej Matlaba

open_system('uklad_IIrz') % Otwarcie modelu Simulinka

wn = 2; % Wartość częstotliwości drgań własnych

zeta = [ 0.25 0.5 0.75]; % Wartości współczynników tłumienia

for i=1:3, % Pętla for

set_param('uklad_IIrz/Transfer Fcn','Numerator',num2str(wn^2),...

'Denominator','[1 2*zeta(i)*wn wn^2]') % Ustawienie

% parametrów transmitancji

sim('uklad_IIrz') % Wykonanie symulacji

t(:,i)=wyniki(:,1); % Podstawienie wyników symulacji

u(:,i)=wyniki(:,2); % pod odpowiednią zmienną

y(:,i)=wyniki(:,3);

end; % Zakończenie pętli for

close_system % Zamknięcie modelu Simulinka

id1 = figure(1) % Otwarcie okna graficznego

plot(t(:,1),u(:,1),'k:',t(:,1),y(:,1),'k-',t(:,2),y(:,2),'k-',...

t(:,3),y(:,3),'k-') % Wykreślenie zmiennych na wykresie

xlabel('t [s]') % Opis osi x

ylabel('y(t)') % Opis osi y

title('Odpowiedzi skokowe układu II rzędu') % Tytuł wykresu

set( id1, 'Color',[1 1 1]) % Usunięcie marginesu otaczającego wykresy

Przykład 4

Dla układu opisanego w przykładzie 1, zamaskować blok z transmitancją II rzędu

]

i ustawić wartości parametrów na ω_n = 1 oraz ζ = 0.5.

Należy zbudować model pokazany na poprzednim rysunku.

Rys. Schemat modelu z zastąpionymi współczynnikami transmitancji na zmienne

Kolejka

Rozważmy kolejkę - zabawkę składającą się z lokomotywy i dołączanego do niej (poprzez sprężysty zaczep) wagonika.

M1 = 1 kg

M2 = 0.5 kg

k = 1 N/sec stała sprężyny

F= 1 N

= 0.002 sec/m. współczynnik tarcia

g = 9.8 m/s^2 przysp. ziemskie

Z praw Newtona wynika, że suma sił działających na ciało jest równa jego masie pomnożonej przez przyspieszenie. W tym wypadku, na ciało M₁działają następujące siły: siła pochodząca od silnika elektrycznego, tarcie oraz siła, którą oddziałuje sprężyna. Na M₂ - jedynie siła sprężyny oraz tarcie. Można to opisać równaniami:

SUWNICA:

Oznaczenia:

• przyspieszenie ziemskie
  

• masę wózka
  

• masę kuli
  

• długość liny
  

• x - położenie środka wózka

• ၪ - wychylenie kuli z położenia równowagi

• F - siła napędowa wytwarzana przez silnik (wymuszenie) - prostokątny impuls o amplitudzie wynoszącej
  i czasie trwania
  

• przyjąć zerowe warunki początkowe (zerowe położenie początkowe i wychylenie jak również prędkości liniowa i kątowa)

 

Równania opisujące model.

 

 

 

 Równania opisujące model są równaniami nieliniowymi (pierwsza pochodna ၪ jest podnoszona do kwadratu, ၪ jest również argumentem funkcji sin i cos). Dlatego układu tego nie można bezpośrednio rozwiązać stosując metodę zmiennych stanu, nie można również wyznaczyć transmitancji. Można go natomiast łatwo rozwiązać tworząc odpowiedni model w Simulinku

 

Układ ten można jednak poddać pewnemu uproszczeniu (linearyzacji). W tym celu utwórzmy sobie następujące zmienne stanu:

 

 Należy przyjąć, że układ pracuje przy małych wychyleniach obciążenia oraz małych prędkościach kątowych tego obciążenia:

cos x₃ Ⴛ 1, sin x₃ Ⴛ x₃, sin² x₃ Ⴛ 0, x₄² Ⴛ 0

wówczas:

 

gdzie:
,
,
,

Modelowanie zawieszenia samochodu:

Równania ruchu:

Poddając stronami powyższe równania transformacie Laplace'a otrzymamy:

 

 

 

 

Po znalezieniu macierzy odwrotnej otrzymamy następujące równanie:

 

 

Kiedy przyjmiemy, że sygnał sterujący U(s) = 0 otrzymamy transmitancję następującej postaci:

---