pobieranie: Modelowanie i Symulacja, ANALIZA WŁASNOŚCI DYNAMICZNYCH WYBRANEGO OBIEKTU FIZYCZNEGO
Pobierz dokument pdfPrzeglądaj wersję html pliku:
Laboratorium Modelowania i symulacji
2008 r.
Wydział Elektryczny Zespół Automatyki (ZTMAiPC) ZERiA
LABORATORIUM MODELOWANIA I SYMULACJI Ćwiczenie 5 ANALIZA WŁASNOŚCI DYNAMICZNYCH WYBRANEGO OBIEKTU FIZYCZNEGO
1. Opis właściwości dynamicznych obiektu
Typowym obiektem dynamicznym, którego opis moŜna przedstawić za pomocą liniowego równania róŜniczkowego jest układ zawieszenia pojazdu. Model blokowy tego układu przedstawia rysunek 1. Pojazd o masie m jest zawieszony nad profilem drogi za pośrednictwem dwóch elementów posiadających właściwości dynamiczne. Pierwszym z nich jest spręŜyna, dla której zaleŜność między siłą do niej przyłoŜoną, a jej odkształceniem wyraŜa się stałym współczynnikiem k. Drugim elementem jest tłumik olejowy, dla którego siła oporu przemieszczenia tłoka i prędkość jego przemieszczania powiązane są stałym współczynnikiem B.
pojazd
m Y=Y0+y(t) k B
profil drogi
a
poziom odniesienia
U=U0+u(t)
Rys.1 . Model blokowy mechanicznego układu zawieszenia samochodu m - masa pojazdu, k – współczynnik spręŜystości zawieszenia, B – prędkościowy współczynnik tłumienia. Po zsumowaniu wszystkich sił działających w układzie i przyjęciu poziomu odniesienia model układu z rys.1 moŜna opisać równaniem:
m
d Y t dt 2
2
B
dY t dt
dU t dt
k Y t
U t
mg ,
(1)
W stanie spoczynku połoŜenie drogi nie zmienia się (przyjmijmy je na poziomie U0). Nie zmienia się równieŜ połoŜenie nadwozia (przyjmijmy je na poziomie Y0). W stanie spoczynku takŜe pochodne są równe zero, dlatego moŜna zapisać:
k Y0 U0
mg ,
(2)
Podstawiając Y=Y0+y(t), U=U0+u(t) oraz wzór (2) do równania (1) otrzymamy zaleŜność zmian połoŜenia nadwozia od zmian profilu drogi:
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-1-
Laboratorium Modelowania i symulacji
2008 r.
m
d y t dt 2
2
B
dy t dt
k y t
B
du t dt
k u t ,
(3)
Równie (3) opisuje dynamikę zawieszenia pojazdu.
2. Cel analizy modelu
Zadaniem zawieszenia jest kompensacja zmian profilu drogi u(t) tak, aby odczuwalne przemieszczenie nadwozia pojazdu y(t) dla pasaŜerów było moŜliwie jak najmniejsze. W celu uproszenia analizy matematycznej przyjmiemy jeden z gorszych przypadków, kiedy droga zmienia profil skokowo o wartość a (Rys. 1). Wówczas zadaniem projektowym jest taki dobór współczynników k oraz B aby przebiegi czasowe odpowiedzi zawieszenia na wymuszenie skokowe o wysokości a były w pewien sposób optymalne. W tym przypadku jako kryterium optymalne przyjmiemy brak oscylacji nadwozia.
3. Program ćwiczenia 3.1. Opis układu w postaci transmitancji operatorowej
1. Na podstawie równania (3) wyznaczyć transmitancję operatorową układu G(s)=Y(s)/X(s). 2. Zasymulować działanie układu zawieszenia wprowadzając skrypt zawieszenie1.m:
clear all clc m=1000 B=200 k=20000 a=0.05 licznik=[B/k,1]; mianownik=[m/k, B/k, 1] sys=tf(licznik,mianownik) figure(10) pzmap(sys) hold on figure(11) step(a*sys) % % % % mas pojazdu[kg] współczynnik tłumienia [Ns/m] współczynnik spręŜystości[N/m] skokowa zmiana profilu drogi o 5 cm
% definicja systemu % rozlokowanie zer i biegunów
% odpowiedz skokowa
3. Czy odpowiedź układu ma charakter oscylacyjny? 4. Jaki wcześniej poznany system dynamiczny ma podobną odpowiedź skokową?
3.2. Warunki oscylacji
Charakter odpowiedzi skokowej systemu zaleŜy od pierwiastków równania charakterystycznego, czyli w naszym przypadku od biegunów transmitancji układu. 1. Przyjmując m=1000 kg i k=20000 N/m wyznaczyć dla jakiej wartości B bieguny transmitancji badanego układu przyjmą wartości rzeczywiste. 2. Zasymulować działanie układu zawieszenia wprowadzając skrypt zawieszenie2.m:
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-2-
Laboratorium Modelowania i symulacji clear all clc m=1000 k=20000 a=0.05 B=sqrt(4*k*m) licznik=[B/k,1]; mianownik=[m/k, B/k, 1] sys=tf(licznik,mianownik) %definicja systemu figure(10) pzmap(sys) % rozlokowanie zer i biegunów hold on figure(11) step(a*sys) %masa pojazdu[kg] %współczynnik spręŜystości[N/m] % skokowa zmiana profilu drogi o 5 cm %wartość graniczna B dla warunku oscylacji
2008 r.
% odpowiedz skokowa
3. Przyjmując B∈<3000;10000> wyznaczyć bieguny odpowiedzi skokowej układu. Zasymulować działanie układu zawieszenia wprowadzając skrypt zawieszenie3.m:
set(0,'DefaultAxesFontSize',12) clear all clc figure(10); clf m=1000 %masa pojazdu[kg] k=20000 %współczynnik spręŜystości[N/m] a=0.05 % skokowa zmiana profilu drogi o 5 cm B=3000:3000:10000; for i=1:max(size(B)) disp(['B= ',num2str(B(i))]) sys=tf([B(i)/k,1],[m/k, B(i)/k, 1]) sys_zpk=zpk(sys) [zera,bieguny,wzmocnienie]=zpkdata(sys_zpk,'v') figure(10); pzmap(sys); xlabel(['współczynnik tłumienia B=3000-10000']); gtext(['B=',num2str(B(i))]); hold on; pause(0.2) figure(10+i) step(sys) title(['współczynnik tłumienia B= ',num2str(B(i))]); end figure(10); sgrid;
Na mapie zer i biegunów naleŜy kliknąć obok pojawiających się nowych zer (o). Na charakterystykach skokowych zaznaczyć czas regulacji, przeregulowanie, czas narastania.
3.3. Aproksymacja układem oscylacyjnym
Z uwagi na podobieństwo postaci mianownika transmitancji rozpatrywanego układu z układem oscylacyjnym o transmitancji:
G osc s
2 n
s2 2
n
s
2 n
(4)
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-3-
Laboratorium Modelowania i symulacji
2008 r.
gdzie ζ - względny współczynnik tłumienia, ωn – pulsacja drgań niegasnących [rad/s] (ωn=2πf), moŜna porównać parametry obydwu układów, a następnie wyznaczyć B oraz k na podstawie zadanej wartości ζ i ωn. Po przekształceniu transmitancji badanego systemu, tak aby wyraz wolny przy najwyŜszej potędze mianownika był równy 1 otrzymamy zaleŜności:
2 ⋅ϖ n ⋅ ξ =
B k , ϖn = m m
(5)
Podejście takie, mimo Ŝe nie jest precyzyjne pozwala na oszacowanie parametrów B, k układu zawieszenia w sposób przybliŜony na podstawie łatwo interpretowalnych parametrów ζ i ωn. 1. Porównać odpowiedzi skokowe, rozkład zer i biegunów układu zawieszenia i układu oscylacyjnego wprowadzając skrypt zawieszenie4.m:
clear all clc m=1000 k=20000 zeta=.2:.6:2; omegan=sqrt(k/m) for i=1:max(size(zeta)) B=2*m*zeta(i)*omegan; disp(['zeta= ',num2str(zeta(i))]) disp(['B= ',num2str(B)]) % aproksymacja parametrów układu zawieszenia układem oscylacyjnym disp('Układ zawieszenia ') sys=tf([B/k,1],[m/k, B/k, 1]) sys_zpk=zpk(sys) [zera,bieguny,wzmocnienie]=zpkdata(sys_zpk,'v') figure(10); pzmap(sys); xlabel(['współczynnik tłumienia (zawieszenie) \zeta= ',num2str(zeta(i))]); gtext(['\zeta=',num2str(zeta(i))]); hold on; pause(0.2) figure(10+i) step(sys) title(['współczynnik tłumienia (zawieszenie) \zeta= ',num2str(zeta(i))]); %układ oscylacyjny disp('Układ oscylacyjny') sys2=tf([omegan],[1, 2*zeta(i)*omegan, omegan*omegan]) sys2_zpk=zpk(sys2) [zera2,bieguny2,wzmocnienie2]=zpkdata(sys2_zpk,'v') figure(40); pzmap(sys2); xlabel(['współczynnik tłumienia (układ oscylacyjny)\zeta= ',num2str(zeta(i))]); gtext(['\zeta=',num2str(zeta(i))]); hold on; pause(0.2) figure(40+i) step(sys2) title(['współczynnik tłumienia (układ oscylacyjny) \zeta= ',num2str(zeta(i))]); end figure(10); sgrid; figure(40); sgrid;
Na mapie zer i biegunów naleŜy kliknąć obok pojawiających się nowych zer (o) lub biegunów(x).
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-4-
Laboratorium Modelowania i symulacji
2008 r.
2. Dla jakich wartości ζ róŜnice w odpowiedziach czasowych są znaczące? (na wykresach zaznaczyć czas regulacji, przeregulowanie, czas narastania) 3. Ocenić połoŜenie zer i biegunów dla transmitancji układu zawieszenia w funkcji ζ. 4. Jaki efekt wywołuje zbliŜenie zera z biegunem?
3.4. Charakterystyki częstotliwościowe układu
1. Porównać charakterystyki częstotliwościowe układu zawieszenia i układu oscylacyjnego wprowadzając skrypt zawieszenie5.m:
clear all clc m=1000; k=20000; zeta=.2:0.6:2; czas=0:.1:2; for i=1:max(size(zeta)) omegan=sqrt(k/m); B=2*m*zeta(i)*omegan; % aproksymacja parametrów układu zawieszenia układem oscylacyjnym sys=tf([B/k,1],[m/k, B/k, 1]); figure(10); pzmap(sys); xlabel(['współczynnik tłumienia (zawieszenie) \zeta= ',num2str(zeta(i))]); gtext(['\zeta=',num2str(zeta(i))]); hold on; pause(0.2) figure(10+i) grid on bode(sys) grid on title(['współczynnik tłumienia (zawieszenie) \zeta= ',num2str(zeta(i))]); %układ oscylacyjny sys2=tf([omegan],[1, 2*zeta(i)*omegan, omegan*omegan]); figure(40); pzmap(sys2); xlabel(['współczynnik tłumienia (układ oscylacyjny)\zeta= ',num2str(zeta(i))]); gtext(['\zeta=',num2str(zeta(i))]); hold on; pause(0.2) figure(40+i) grid on bode(sys2) grid on title(['współczynnik tłumienia (układ oscylacyjny) \zeta= ',num2str(zeta(i))]); end figure(10); sgrid; figure(40); sgrid;
Na mapie zer i biegunów naleŜy kliknąć obok pojawiających się nowych zer (o) lub biegunów(x).
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-5-
Modelowanie i Symulacja, ANALIZA WŁASNOŚCI DYNAMICZNYCH WYBRANEGO OBIEKTU FIZYCZNEGO
Laboratorium Modelowania i symulacji
2008 r.
Wydział Elektryczny Zespół Automatyki (ZTMAiPC) ZERiA
LABORATORIUM MODELOWANIA I SYMULACJI Ćwiczenie 5 ANALIZA WŁASNOŚCI DYNAMICZNYCH WYBRANEGO OBIEKTU FIZYCZNEGO
1. Opis właściwości dynamicznych obiektu
Typowym obiektem dynamicznym, którego opis moŜna przedstawić za pomocą liniowego równania róŜniczkowego jest układ zawieszenia pojazdu. Model blokowy tego układu przedstawia rysunek 1. Pojazd o masie m jest zawieszony nad profilem drogi za pośrednictwem dwóch elementów posiadających właściwości dynamiczne. Pierwszym z nich jest spręŜyna, dla której zaleŜność między siłą do niej przyłoŜoną, a jej odkształceniem wyraŜa się stałym współczynnikiem k. Drugim elementem jest tłumik olejowy, dla którego siła oporu przemieszczenia tłoka i prędkość jego przemieszczania powiązane są stałym współczynnikiem B.
pojazd
m Y=Y0+y(t) k B
profil drogi
a
poziom odniesienia
U=U0+u(t)
Rys.1 . Model blokowy mechanicznego układu zawieszenia samochodu m - masa pojazdu, k – współczynnik spręŜystości zawieszenia, B – prędkościowy współczynnik tłumienia. Po zsumowaniu wszystkich sił działających w układzie i przyjęciu poziomu odniesienia model układu z rys.1 moŜna opisać równaniem:
m
d Y t dt 2
2
B
dY t dt
dU t dt
k Y t
U t
mg ,
(1)
W stanie spoczynku połoŜenie drogi nie zmienia się (przyjmijmy je na poziomie U0). Nie zmienia się równieŜ połoŜenie nadwozia (przyjmijmy je na poziomie Y0). W stanie spoczynku takŜe pochodne są równe zero, dlatego moŜna zapisać:
k Y0 U0
mg ,
(2)
Podstawiając Y=Y0+y(t), U=U0+u(t) oraz wzór (2) do równania (1) otrzymamy zaleŜność zmian połoŜenia nadwozia od zmian profilu drogi:
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-1-
Laboratorium Modelowania i symulacji
2008 r.
m
d y t dt 2
2
B
dy t dt
k y t
B
du t dt
k u t ,
(3)
Równie (3) opisuje dynamikę zawieszenia pojazdu.
2. Cel analizy modelu
Zadaniem zawieszenia jest kompensacja zmian profilu drogi u(t) tak, aby odczuwalne przemieszczenie nadwozia pojazdu y(t) dla pasaŜerów było moŜliwie jak najmniejsze. W celu uproszenia analizy matematycznej przyjmiemy jeden z gorszych przypadków, kiedy droga zmienia profil skokowo o wartość a (Rys. 1). Wówczas zadaniem projektowym jest taki dobór współczynników k oraz B aby przebiegi czasowe odpowiedzi zawieszenia na wymuszenie skokowe o wysokości a były w pewien sposób optymalne. W tym przypadku jako kryterium optymalne przyjmiemy brak oscylacji nadwozia.
3. Program ćwiczenia 3.1. Opis układu w postaci transmitancji operatorowej
1. Na podstawie równania (3) wyznaczyć transmitancję operatorową układu G(s)=Y(s)/X(s). 2. Zasymulować działanie układu zawieszenia wprowadzając skrypt zawieszenie1.m:
clear all clc m=1000 B=200 k=20000 a=0.05 licznik=[B/k,1]; mianownik=[m/k, B/k, 1] sys=tf(licznik,mianownik) figure(10) pzmap(sys) hold on figure(11) step(a*sys) % % % % mas pojazdu[kg] współczynnik tłumienia [Ns/m] współczynnik spręŜystości[N/m] skokowa zmiana profilu drogi o 5 cm
% definicja systemu % rozlokowanie zer i biegunów
% odpowiedz skokowa
3. Czy odpowiedź układu ma charakter oscylacyjny? 4. Jaki wcześniej poznany system dynamiczny ma podobną odpowiedź skokową?
3.2. Warunki oscylacji
Charakter odpowiedzi skokowej systemu zaleŜy od pierwiastków równania charakterystycznego, czyli w naszym przypadku od biegunów transmitancji układu. 1. Przyjmując m=1000 kg i k=20000 N/m wyznaczyć dla jakiej wartości B bieguny transmitancji badanego układu przyjmą wartości rzeczywiste. 2. Zasymulować działanie układu zawieszenia wprowadzając skrypt zawieszenie2.m:
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-2-
Laboratorium Modelowania i symulacji clear all clc m=1000 k=20000 a=0.05 B=sqrt(4*k*m) licznik=[B/k,1]; mianownik=[m/k, B/k, 1] sys=tf(licznik,mianownik) %definicja systemu figure(10) pzmap(sys) % rozlokowanie zer i biegunów hold on figure(11) step(a*sys) %masa pojazdu[kg] %współczynnik spręŜystości[N/m] % skokowa zmiana profilu drogi o 5 cm %wartość graniczna B dla warunku oscylacji
2008 r.
% odpowiedz skokowa
3. Przyjmując B∈<3000;10000> wyznaczyć bieguny odpowiedzi skokowej układu. Zasymulować działanie układu zawieszenia wprowadzając skrypt zawieszenie3.m:
set(0,'DefaultAxesFontSize',12) clear all clc figure(10); clf m=1000 %masa pojazdu[kg] k=20000 %współczynnik spręŜystości[N/m] a=0.05 % skokowa zmiana profilu drogi o 5 cm B=3000:3000:10000; for i=1:max(size(B)) disp(['B= ',num2str(B(i))]) sys=tf([B(i)/k,1],[m/k, B(i)/k, 1]) sys_zpk=zpk(sys) [zera,bieguny,wzmocnienie]=zpkdata(sys_zpk,'v') figure(10); pzmap(sys); xlabel(['współczynnik tłumienia B=3000-10000']); gtext(['B=',num2str(B(i))]); hold on; pause(0.2) figure(10+i) step(sys) title(['współczynnik tłumienia B= ',num2str(B(i))]); end figure(10); sgrid;
Na mapie zer i biegunów naleŜy kliknąć obok pojawiających się nowych zer (o). Na charakterystykach skokowych zaznaczyć czas regulacji, przeregulowanie, czas narastania.
3.3. Aproksymacja układem oscylacyjnym
Z uwagi na podobieństwo postaci mianownika transmitancji rozpatrywanego układu z układem oscylacyjnym o transmitancji:
G osc s
2 n
s2 2
n
s
2 n
(4)
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-3-
Laboratorium Modelowania i symulacji
2008 r.
gdzie ζ - względny współczynnik tłumienia, ωn – pulsacja drgań niegasnących [rad/s] (ωn=2πf), moŜna porównać parametry obydwu układów, a następnie wyznaczyć B oraz k na podstawie zadanej wartości ζ i ωn. Po przekształceniu transmitancji badanego systemu, tak aby wyraz wolny przy najwyŜszej potędze mianownika był równy 1 otrzymamy zaleŜności:
2 ⋅ϖ n ⋅ ξ =
B k , ϖn = m m
(5)
Podejście takie, mimo Ŝe nie jest precyzyjne pozwala na oszacowanie parametrów B, k układu zawieszenia w sposób przybliŜony na podstawie łatwo interpretowalnych parametrów ζ i ωn. 1. Porównać odpowiedzi skokowe, rozkład zer i biegunów układu zawieszenia i układu oscylacyjnego wprowadzając skrypt zawieszenie4.m:
clear all clc m=1000 k=20000 zeta=.2:.6:2; omegan=sqrt(k/m) for i=1:max(size(zeta)) B=2*m*zeta(i)*omegan; disp(['zeta= ',num2str(zeta(i))]) disp(['B= ',num2str(B)]) % aproksymacja parametrów układu zawieszenia układem oscylacyjnym disp('Układ zawieszenia ') sys=tf([B/k,1],[m/k, B/k, 1]) sys_zpk=zpk(sys) [zera,bieguny,wzmocnienie]=zpkdata(sys_zpk,'v') figure(10); pzmap(sys); xlabel(['współczynnik tłumienia (zawieszenie) \zeta= ',num2str(zeta(i))]); gtext(['\zeta=',num2str(zeta(i))]); hold on; pause(0.2) figure(10+i) step(sys) title(['współczynnik tłumienia (zawieszenie) \zeta= ',num2str(zeta(i))]); %układ oscylacyjny disp('Układ oscylacyjny') sys2=tf([omegan],[1, 2*zeta(i)*omegan, omegan*omegan]) sys2_zpk=zpk(sys2) [zera2,bieguny2,wzmocnienie2]=zpkdata(sys2_zpk,'v') figure(40); pzmap(sys2); xlabel(['współczynnik tłumienia (układ oscylacyjny)\zeta= ',num2str(zeta(i))]); gtext(['\zeta=',num2str(zeta(i))]); hold on; pause(0.2) figure(40+i) step(sys2) title(['współczynnik tłumienia (układ oscylacyjny) \zeta= ',num2str(zeta(i))]); end figure(10); sgrid; figure(40); sgrid;
Na mapie zer i biegunów naleŜy kliknąć obok pojawiających się nowych zer (o) lub biegunów(x).
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-4-
Laboratorium Modelowania i symulacji
2008 r.
2. Dla jakich wartości ζ róŜnice w odpowiedziach czasowych są znaczące? (na wykresach zaznaczyć czas regulacji, przeregulowanie, czas narastania) 3. Ocenić połoŜenie zer i biegunów dla transmitancji układu zawieszenia w funkcji ζ. 4. Jaki efekt wywołuje zbliŜenie zera z biegunem?
3.4. Charakterystyki częstotliwościowe układu
1. Porównać charakterystyki częstotliwościowe układu zawieszenia i układu oscylacyjnego wprowadzając skrypt zawieszenie5.m:
clear all clc m=1000; k=20000; zeta=.2:0.6:2; czas=0:.1:2; for i=1:max(size(zeta)) omegan=sqrt(k/m); B=2*m*zeta(i)*omegan; % aproksymacja parametrów układu zawieszenia układem oscylacyjnym sys=tf([B/k,1],[m/k, B/k, 1]); figure(10); pzmap(sys); xlabel(['współczynnik tłumienia (zawieszenie) \zeta= ',num2str(zeta(i))]); gtext(['\zeta=',num2str(zeta(i))]); hold on; pause(0.2) figure(10+i) grid on bode(sys) grid on title(['współczynnik tłumienia (zawieszenie) \zeta= ',num2str(zeta(i))]); %układ oscylacyjny sys2=tf([omegan],[1, 2*zeta(i)*omegan, omegan*omegan]); figure(40); pzmap(sys2); xlabel(['współczynnik tłumienia (układ oscylacyjny)\zeta= ',num2str(zeta(i))]); gtext(['\zeta=',num2str(zeta(i))]); hold on; pause(0.2) figure(40+i) grid on bode(sys2) grid on title(['współczynnik tłumienia (układ oscylacyjny) \zeta= ',num2str(zeta(i))]); end figure(10); sgrid; figure(40); sgrid;
Na mapie zer i biegunów naleŜy kliknąć obok pojawiających się nowych zer (o) lub biegunów(x).
Ćwiczenie 5 – Analiza własności dynamicznych wybranego obiektu fizycznego
-5-
