Zagadnienia

6.1 Bezstresowe rozwiązanie problemu
- 6.1.1 Próba weryfikacji rozwiązania
6.2 Czy interpretacja wyniku jest poprawna?
- 6.2.1 Jak dobrać dobrą rozdzielczość?
6.3 Eksperymenty z parametrami układu
6.4 Uwagi i uzupełnienia

6. Charakterystyka pracy transformatora

Do transformatora o $N$ zwojach doprowadzony jest prąd zmienny o napięciu $E$ i częstotliwości $\omega/2\pi$ . W modelowaniu jego pracy [7, Chapter 6, APP10] wykorzystuje się równanie różniczkowe opisujące zmiany natężenia prądu $\phi$ w czasie:

$\phi^{{\prime\prime}}+\omega _{0}^{2}\phi+b\phi^{3}=\frac{\omega}{N}E\cos(\omega t).$

(6.1)

Parametry $\omega _{0}$ oraz $b$ zależą od konkretnego transformatora. Podobne równanie rozwiązuje się, modelując np. zjawisko tzw. ferrorezonansu w dużych sieciach transformatorów wysokiego napięcia, lecz wtedy czasem człon nieliniowy $\phi^{3}$ zastępuje się wyższą nieparzystą potęgą $\phi$ , np. $\phi^{{11}}$ [18].

Niech warunek początkowy będzie postaci $\phi(0)=\phi^{{\prime}}(0)=0$ . Naszym celem jest

narysowanie wykresu $\phi(t)$ dla $t\in[0,5]$ dla przypadku $E=165$ , $\omega=120\pi$ , $N=600$ , $\omega _{0}^{2}=83$ , $b=0.14$ ;
zbadanie wpływu niewielkich zmian parametru $b$ na przebieg rozwiązania.

6.1. Bezstresowe rozwiązanie problemu

Aby rozwiązać (6.1), zastosujemy standardowy solver równań różniczkowych z Octave, lsode. Po sprowadzeniu równania do układu dwóch równań pierwszego rzędu,

$\dfrac{d}{dt}\begin{pmatrix}\phi\\ v\end{pmatrix}=\begin{pmatrix}v\\ -\omega _{0}^{2}\phi-b\phi^{3}+\dfrac{\omega}{N}E\cos(\omega t)\end{pmatrix},$

określamy funkcję wyznaczającą wektor prawej strony (przy umowie, że $y=(y_{1},y_{2})=(\phi,v)$ ):

E = 165; omega = 120*pi; N = 600; omega02 = 83; b = 0.14;
oEN = (omega*E)/N;

transf = @(y,t) [y(2); -y(1)*(omega02 + b*y(1)^2) + oEN*cos(omega*t)];

Numeryczne rozwiązanie i wykres $\phi$ (czyli pierwszej współrzędnej wektora $y$ ) wyznaczymy stosując proste komendy:

y0 = [0;0]; % warunek początkowy

K = 100; % rozdzielczość wykresu
t = linspace(0,5,K);

y = lsode(transf, y0, t);
plot(t, y(:,1)); % rysujemy tylko pierwszą współrzedną

Ładny, regularny i periodyczny wykres sugeruje, że --- na mocy ,,estetycznego kryterium weryfikacji wyników'' --- rozwiązanie jest poprawne i sensowne.

Rys. 6.1. Rozwiązanie równania transformatora? To się jeszcze okaże.

6.1.1. Próba weryfikacji rozwiązania

Uzyskany przez nas na rysunku 6.1 wykres — regularny, ,,optycznie prawie periodyczny” — zdaje się sensowny. Jednak — pamiętając o tym, że niektóre zadania mogą być numerycznie trudne, a ich rozwiązania zwodniczo eleganckie — na wszelki wypadek zastosujemy inżynierski chwyt, polegający na weryfikacji rozwiązania poprzez drastyczne zaostrzenie kryteriów tolerancji błędu.

Aby to osiągnąć, musimy nieco pomanipulować przy parametrach pracy lsode. Najpierw zapamiętamy bieżące wartości tolerancji błędu¹⁴Przypomnijmy na wszelki wypadek, że ustalenie wartości tych parametrów nie spowoduje bynajmniej, że uzyskane rozwiązanie będzie obarczone właśnie takim błędem; faktyczna rola parametrów tolerancji błędu jest opisana specjalistycznej literaturze, np. w [9]. względnego i bezwzględnego:

rtol = lsode_options('relative tolerance');
atol = lsode_options('absolute tolerance');

Następnie ustalimy nowe wartości na poziomie $10^{4}$ razy(!) mniejszym:

lsode_options('relative tolerance', rtol*1e-4);
lsode_options('absolute tolerance', atol*1e-4);

i na nowo uruchomimy tą samą symulację:

y = lsode(transf, y0, t);
plot(t, y(:,1));

% na zakończenie przywracamy poprzednie wartości parametrów lsode
lsode_options('relative tolerance', rtol);
lsode_options('absolute tolerance', atol);

Choć na wyniki przyjdzie nam czekać znacznie dłużej niż poprzednio — wszak utrudniliśmy pracę solverowi! — efekt, który możemy sprawdzić na rysunku 6.2, dość przekonująco potwierdza, że solver nie oszukał nas poprzednim razem i wyznaczone rozwiązanie jest praktycznie identyczne.

$Po drastycznym zredukowaniu tolerancji błędu lsode, wykres praktycznie nie zmienia się --- to znak, że \en{solver} prawdopodobnie działa dobrze.$

Rys. 6.2. Numeryczne potwierdzenie jakości rozwiązania równania transformatora? ,,Na oko” zdaje się, że jest OK.

Zatem na podstawie uzyskanego wykresu możemy skonkludować, że nasz transformator ma na tyle dużą bezwładność, że — choć w jego pracy pojawiają się oscylacje — przebiegają one łagodnie. W badanym okresie czasu natężenie tylko dwa razy (pięć, jeśli darować sobie drobne różnice) przyjęło maksymalną wartość.

6.2. Czy interpretacja wyniku jest poprawna?

Jednak otrzymany przez nas wykres jest — w potocznym rozumieniu — całkowicie błędny! Zapomnieliśmy bowiem o innym zjawisku: o aliasingu, czyli niebezpieczeństwie wyciągnięcia błędnych wniosków o przebiegu funkcji w przypadku zbyt małej liczby punktów próbkowania jej wartości. Rzeczywiście, powinni byliśmy na wszelki wypadek zbadać, czy zwiększenie rozdzielczości wykresu nie spowoduje jego wyraźnej zmiany….

Zatem zobaczmy. Wybierzemy — zamiast $K=100$ punktów wykresu, jak dotychczas — na początek $K=300$ równoodległych punktów:

K = 300;
t = linspace(0,5,K);

y = lsode(transf, y0, t);
plot(t, y(:,1));

Wykres jest zupełnie inny niż dla K=100.

Rys. 6.3. Rozwiązanie równania transformatora próbkowane w $K=300$ punktach.

Hmmm… dostaliśmy wykres (zob. rysunek 6.3), w żaden sposób nie przypominający rysunku 6.1, który jeszcze przed chwilą wydawał nam się tak sensowny… A przecież jedynym zmienionym parametrem była liczba punktów próbkowania rozwiązania: zatem diagnoza jest jasna: faktycznie, mamy aliasing.

No dobrze, zatem który z wykresów jest prawidłowy? Możemy powtórzyć nasz chwyt i próbkować w $K=600$ punktach (zob. rysunek 6.4), ale… może lepiej postąpić nieco bardziej metodycznie, zwłaszcza, że za chwilę możemy natknąć się na barierę rozdzielczości… naszego ekranu komputerowego?

Wykres jest zupełnie inny niż dla K=100 i innny niż dla K=300...

Rys. 6.4. Rozwiązanie równania transformatora próbkowane w $K=600$ punktach.

Już przekonaliśmy się, że (zapewne) nasz solver, lsode, nie ma większych kłopotów z aproksymacją rozwiązania. Pozostaje więc przekonać się, kiedy niewielkie zwiększanie (lub, co — w świetle powyższego — na jedno wychodzi, zmniejszenie) rozdzielczości wykresu przestanie powodować istotną zmianę wykresu rozwiązania.

6.2.1. Jak dobrać dobrą rozdzielczość?

Gdybyśmy bowiem wcześniej choć chwilę zastanowili się nad charakterem naszego równania, zobaczylibyśmy, że jego prawa strona — zewnętrzne wymuszenie — jest postaci $\cos(\omega\, t)$ , czyli, dla naszych danych, postaci $\cos(120\,\pi\, t)$ . Nawet intuicyjnie jest jasne, że próbkowanie wartości rozwiązania powinniśmy dopasować do tej częstości. W interesującym nas przedziale czasowym $t\in[0,5]$ takie wymuszenie będzie miało $120\cdot 5=700$ ekstremów, więc — biorąc skromnie po 10 węzłów pomiędzy jednym a drugim ekstremum — powinniśmy użyć około 7000 węzłów, by ustawić wstępną (minimalną sensowną?) rozdzielczość wizualizacji.

Naturalnie, możemy dodatkowo oprzeć się na naszej wiedzy dla liniowych równań różniczkowych. W przypadku $b=0$ drgania będą miały dwie składowe: wolnozmienną postaci $\cos(\omega _{0}t)$ (być może z jakimś przesunięciem fazowym), związaną z rozwiązaniem ogólnym jednorodnego równania oscylatora harmonicznego — oraz składową szybkozmienną, postaci $\cos(\omega t)$ , będącą rozwiązaniem szczególnym równania niejednorodnego. W naszym przypadku, da to coś w stylu $A_{0}\cos(9.11\, t+\phi)+A_{1}\cos(376.99\, t)$ — zatem faktycznie musimy dostosować rozdzielczość wizualizacji co najmniej do najszybciej zmieniającego się składnika, $\cos(\omega t)$ .

Aby utwierdzić się w naszych przypuszczeniach, zrealizujemy w pętli wyświetlanie rozwiązania, oryginalnie wypróbkowanego w bardzo wielu punktach ( $K=5000$ ), nie zważając na rozdzielczość monitora, która zapewne nie przekracza 2000 pikseli w poziomie, na coraz rzadszej siatce węzłów:

K = 5000;
t = linspace(0,T,K); % highest possible resolution
y = lsode(transf, y0, t); % reference solution
plot(t,y(:,1)); title('Reference plot'); pause(3);

P = 1;
for k = [1:10]
	plot(t(P:k:end),y(P:k:end,1));
	title(['Reference plot for K = ',num2str(K),', every ',num2str(k),'th point']); pause(3);
end

Wykres nieco się ustabilizował i pokazuje wyraźne nakładanie się dwóch składowych o różnych częstotliwościach.

Rys. 6.5. Rozwiązanie równania transformatora próbkowane w $K=5000$ punktach.

Ta seria wykresów przekonuje nas jasno, że rozwiązanie zmienia się zbyt szybko, by było czytelnie wizualizowane na odcinku $[0,5]$ — zob. rysunek 6.5. Ponadto faktycznie wydaje się, że obserwujemy zjawisko nakładania się szybkozmiennej składowej na składową wolnozmienną. Obie te obserwacje skłaniają więc nas do tego, by przyjrzeć się dokładniej wykresowi ograniczonemu do znacznie węższego zakresu $t$ : wybieramy ,,na oko” okno o szerokości $5\%$ oryginalnego przedziału, zlokalizowane na jego końcu.

P = floor(0.95*K); % ostatnie 5% wykresu
for k = [1:10]
	plot(t(P:k:end),y(P:k:end,1));
	title(['Reference plot for K = ',num2str(K),', every ',num2str(k),'th point']); pause(3);
end

Wygląda na to, że dysponujemy dostateczną (niewiele większą od minimalnej) liczbą punktów próbkowania w stosunku do realnego przebiegu rozwiązania.

Rys. 6.6. Rozwiązanie równania transformatora próbkowane w $K=5000$ punktach, zobrazowane na krótkim podprzedziale.

Strzał okazał się znakomity (rysunek 6.6) — wydaje się, że przy oryginalnej rozdzielczości $K=5000$ (której w zawężonym przedziale długości $0.25$ odpowiadało 250 węzłów) udało się zlokalizować najważniejsze cechy rozwiązania: wszystko bowiem wskazuje na to, że rozwiązanie ma w badanym zawężonym przedziale około 20 równomiernie rozmieszczonych ekstremów.

Możemy teraz potwierdzić swoje przypuszczenie, biorąc $K=20000$ (czemu odpowiada 1000 węzłów w zawężonym przedziale) i sprawdzając tam przebieg rozwiązania.

Ćwiczenie 6.1

Jak oglądać ostatnie 5% rozwiązania, ale bez $K$ idącego w tysiące (dziesiątki tysięcy)?

E = 165; omega = 120*pi; N = 600; omega02 = 83; b = 0.14;
oEN = (omega*E)/N;

transf = @(y,t) [y(2); -y(1)*(omega02 + b*y(1)^2) + oEN*cos(omega*t)];

y0 = [0;0]; T = 5;

% high voltage transformer
% E = 0.15*63.5e3; C = 777; R = 48.4e3; 
% omegaS = 377/(2*pi); omega02 = (53/(2*pi))^2; omega22 = (85/(2*pi))^2;
% duff = @(y,t) [y(2); ...
% omegaS*E*cos(omegaS*t)-y(1)*omega02-omega22*y(1)^3-y(2)/(R*C) ];

K = 100;
t = linspace(0,T,K);

rtol = lsode_options('relative tolerance');
atol = lsode_options('absolute tolerance');

y = lsode(transf,y0, t);
plot(t,y(:,1)); pause(3);
print('-depsc','duffing-aliasing.eps');

lsode_options('relative tolerance', rtol*1e-4);
lsode_options('absolute tolerance', atol*1e-4);

y = lsode(transf,y0, t); 
plot(t,y(:,1)); pause(3);
print('-depsc','duffing-aliasing-hires.eps');

lsode_options('absolute tolerance', atol);
lsode_options('relative tolerance', rtol);

K = 300;
t = linspace(0,T,K);
y = lsode(transf,y0, t); 
plot(t,y(:,1)); pause(3);
print('-depsc','duffing-aliasing-300.eps');

K = 600;
t = linspace(0,T,K);
y = lsode(transf,y0, t); 
plot(t,y(:,1)); pause(3);
print('-depsc','duffing-aliasing-600.eps');

K = 5000;
t = linspace(0,T,K); % highest possible resolution
y = lsode(transf, y0, t); % reference solution
plot(t,y(:,1)); title('Reference plot'); pause(3);
print('-depsc','duffing-reference.eps');

P = 1;
for k = [1:10]
	plot(t(P:k:end),y(P:k:end,1)); 
	title(['Reference plot for K = ',num2str(K),', every ',num2str(k),'th point']); pause(1);
end

P = floor(0.95*K); % ostatnie 5% wykresu 
for k = [1:10]
	plot(t(P:k:end),y(P:k:end,1)); 
	title(['Reference plot for K = ',num2str(K),', every ',num2str(k),'th point']); pause(1);
end

plot(t(P:end),y(P:end,1));
print('-depsc','duffing-reference-zoom20.eps');

fy = fft(y(:,1)); % fourier transform on phi
fy(1) = []; % pierwsza wspolrzedna to suma wartosci y i jest nieistotna
nyq = T/2;
freq = (1:K/2)/(K/2)*nyq;
semilogx(freq,abs(fy(1:K/2)));

% verifyduffing;

6.3. Eksperymenty z parametrami układu

Naszym drugim celem jest zbadanie wrażliwości rozwiązania na małe zmiany $b$ . Otóż korzystając z poprzednio opracowanego kodu i stopniowo zmniejszając $b$ aż do $0$ widzimy, że rozwiązanie różni się niewiele od dotychczas rozważanego. Możemy więc przypuszczać, że (dla małych wartości $b$ ) rozwiązania zależą w sposób ciągły od tego parametru, a linearyzacja modelu jest tu najprawdopodobniej uzasadniona.

Tu kończy się etap obliczeń naukowych — uprawiania ,,matematyki” eksperymentalnej: teraz powinniśmy przejść do dowodzenia twierdzeń. Czy jesteś w stanie udowodnić powyższe własności rozważanego równania?

6.4. Uwagi i uzupełnienia

Warto na zakończenie wspomnieć, że — choć sprawiły nam niejakie kłopoty — rozwiązania naszego równania przy danych wartościach parametrów zachowywały się bardzo przyzwoicie. Tymczasem dla innych zestawów parametrów może być znacznie gorzej. Równanie tej postaci jest znane pod nazwą równania Duffinga (z periodycznym wymuszeniem). Modeluje ono nie tylko pracę transformatora, ale także działanie pewnych układów mechanicznych. W ogólnym przypadku jest ono postaci

$\phi^{{\prime\prime}}+\epsilon\phi^{{\prime}}\pm\omega _{0}^{2}\phi+\delta\phi^{3}=f(t)$

i jest znanym przykładem zadania różniczkowego, które dla pewnych wartości $\epsilon$ i $\omega _{0}$ może mieć chaotycznie zachowujące się rozwiązania, a to oznacza, że nieunikniony błąd numerycznej aproksymacji szybko spowoduje przeskok na odległą od rzeczywistej trajektorię i całkowitą utratę jakości rozwiązania.

Ćwiczenie 6.2

W modelu sieci energetycznej z transformatorami wysokonapięciowymi pojawia się równanie opisujące wartość natężenia prądu $\phi(t)$ w chwili $t$ , postaci [18]

$\phi^{{\prime\prime}}+\dfrac{1}{RC}\phi^{{\prime}}+\omega _{0}^{2}\phi+\omega _{2}^{2}\phi^{3}=\omega E\cos(\omega t).$

Równanie uzupełnione jest warunkiem początkowym $\phi(0)=...,\phi^{{\prime}}(0)=0$ . Dla wartości parametrów: $R=48.4\text{ [k$\Omega$]}$ , $C=777\text{ [nF]}$ , , $E=0.15\cdot 63.5\text{ [kV]}$ , $\omega=377\text{ [rad/s]}$ , $\omega _{0}=53\text{ [rad/s]}$ , $\omega _{2}=85\text{ [rad/s]}$ , wyznacz przebieg $\phi$ w ciągu pierwszych 5 minut.