Карачун В.В., Мельник В.М., Калініна М.Ф.,

Кривець О.О., Саверченко В.Г.

Національний технічний університет України «КПІ»

ЧАСТОТНІ ХАРАКТЕРИСТИКИ КОНСЕРВАТИВНИХ СИСТЕМ

Консервативні системи, тобто системи, в яких відсутнє розсіювання енергії, або, відповідно, з підживленням енергії в міру її розсіювання, відрізняються тією властивістю, що їх коливання, одного разу почавшись, не затухають. Тому, наприклад, експериментальне ( за допомогою фізичного або математичного експерименту) чисельне визначення їх частотних характеристик наштовхується на істотні труднощі.

В той же час, саме для подібних систем визначення їх частотних характеристик часто є абсолютно необхідними для дослідника, бо саме такі системи особливо чутливо реагують на збурення на резонансних, або близьких до них, частотах.

Можна, звичайно, отримати аналітичний вираз для передаточної функції об'єкта, замінивши в ній на і, тим самим, сформувати вирази для частотних характеристик, але для складних динамічних систем цей шлях надзвичайно громіздкий, а що стосується систем з розподіленими параметрами, то він взагалі практично нереалізуємий внаслідок виникаючих суто математичних труднощів.

Разом з тим, методи чисельного розрахунку перехідних процесів, орієнтовані на реалізацію їх на сучасних комп'ютерах, зараз досить добре відпрацьовані. За необхідності обчислення, наприклад, вільних коливань багатомасової механічної системи із зосередженими параметрами, або об'єкта з розподіленими параметрами, особливих проблем не виникає.

Маючи ж масив ординат вихідного сигналу об'єкта, можна провести його спектральний аналіз і визначити для початку його власні (резонансні) частоти.

Математичною базою спектрального аналізу є перетворення Фур’є

, (1)

де є спектральна характеристика функції .

Її можна навести у у вигляді -

, (2)

де

або, відповідно,

, (3)

де

Для практичної реалізації формул для і в (2) межі інтегрування доводиться обирати від 0 до (будемо називати його часом спостереження процесу). Та обставина, що - скінченне, призводить до того, що визначається не дійсна спектральна характеристика (1), а, так званий, поточний спектр. За він буде наближатись до істинного. Якщо є періодичною функцією, тоді в міру зростання буде наближатися до лінійчатого спектру, де екстремуми будуть мати місце на частотах, що наближаються до власних частот об'єкта (системи), якщо відображає вільний рух. Коли ж об'єкт, що досліджується, знаходиться під дією зовнішніх впливів, тоді в його реакції будуть присутні періодичні компоненти, відповідні гармонічним складовим його входів.

Зупинимося спочатку на визначенні резонансних частот консервативного об'єкту. Якою величиною можна обмежитися? Чим більше, тим краще? Так, звичайно, але тут можливі проблеми з обсягом пам'яті, необхідної для зберігання масиву ординат виходу об'єкта, з іншого боку, - зростає час формування подібного масиву (особливо для складних розподілених об'єктів), а потім час його обробки з обчисленням інтегралів за формулами (2) для досить широкого діапазону частот. Пропонується алгоритм формування на екрані комп'ютера графіка поточного спектра вихідного сигналу об'єкта в міру формування цього сигналу, що позбавляє від необхідності накопичувати (зберігати) цей сигнал. Вважаємо, що він формується або в процесі чисельного інтегрування відповідної системи диференціальних рівнянь, що моделює об'єкт дослідження, або ж подається на комп'ютер через відповідні перетворювачі від датчиків, встановлених на натурному (фізичному) об'єкті.

Обчислення інтегралів за формулою (2) починається з деякого початкового значення (задається користувачем), або навіть з . Для того, щоб користувач отримав уявлення про характер зміни на відрізку часу , на екран можна вивести графік, а потім з інтервалом в кроків чисельного інтегрування (задає користувач) виводиться графік таким чином, щоб черговий графік накладався на попередні і користувач мав можливість спостерігати динаміку зміни в міру зростання часу спостереження. Користувачеві надається можливість в будь-який момент часу звільнити екран, якщо там нагромадилося досить багато кривих, які заважають оцінці останніх на поточний момент часу графіків спектральної щільності. Користувач також може в довільному напрямі змінювати масштаб побудови графіка (в міру зростання t , висота піків на графіку зростає і масштаб доцільно збільшити). Після того як екстремуми, що цікавлять користувача, оформились, з його точки зору, досить чітко, він може зупинити процес і надати графіку остаточного вигляду, скориставшись можливостями, які надаються для цього модулем GromU .

Procedure Specter;

Const Tsu='1-стир., 2-змін.граф., 3-збільш.граф., Esc-вихід';

r:integer=100;

var S, U, V:integer;

у, Wz, Wt, W:real;

А, В:CoefL;

Ts:=string[10];

Tss:=string[50];

Procedure Graphic;

Var z, s:integer;

Begin

For z:=0 to L do

Begin

s:=Y0-round(sqrt(sqr(А[z])+sqr(В[z]))/Dy);

If z=0 then MoveTo(Xu, s)

Else LineTo(Xu+z, s)

End;

End;

Begin

CoefN; Xu:=(GetMaxX-L) div 2; Init; Ft[0]:=Yt[0];

PutA; SetColor(1); Rectangle(Xu, 1, Xu+L, 5);

For S:=1 to L do

Begin

Step(ks); Ft[S]:=Yt[NXs]; U:=Xu+S;

Line(Xu+s,2, Xu+s,4)

End;

Puta;

Our(‘Wmax-макс.значення частоти', Wmax);

Puta; Oui(‘R-інтервал оновлення зобр.', R);

Dw:=wmax/L; PutA; Rectangle(Xu, 1, Xu+L, 5);

Ymax:=0; W:=Dw*ks*Tau;

For z:=0 to L do

Begin

Wz:=z*W; А[z]:=0; B[z]:=0;

For s:=0 to L-1 do

Begin

Wt:=Wz*s; А[z]:=А[z]+Ft[s]*cos(wt);

B[z]:=B[z]+Ft[s]*sin(wt)

End;

Y:=sqrt(sqr(А[z])+sqr(В[z]));

If Y>Ymax then Ymax:=Y;

U:=Xu+z; Line(Xu+z,2, Xu+s,4)

End;

Xmin:=0; Xmax:=Wmax; Ymin:=0;

Ymax:=Ymax*5; X0Y0(false);

PutA; ClearDevice; SystCoor; Graphic;

SetViewPort(0,0, GetMaxX, GetMaxY, true);

U:=0; S:=L; V:=0; Ou(Tsu);

Repeat

Step(ks); Y:=Yt[NXs]; inc(U); inc(s);

If s>30000 then

Begin Inc(V); S:=V End;

For z:=0 to L do

Begin

Wt:=z*Dw*t; А[z]:=А[z]+Y*cos(wt);

B[z]:=B[z]+Y*sin(wt)

End;

If U=r then

Begin

Str(t:10:3, Ts); Tss:=' Т='+Ts;

Str(S, Ts); Tss:=Tss+', S='+Ts; PutA;

OutTextXY(420,1, Tss); Ou(Tsu); Graphic; U:=0;

If KeyPressed then J:=ReadKey;

If J<>#27 then

Begin

Case J of

'2': Ymax:=Ymax*1.2;

'3': Ymax:=Ymax/1.2

end;

if J in ['1'..'3'] then

begin

X0Y0(false); J:=' N'; ClearDevice;

SystCoor; Info; Graphic; Ou(Tsu);

OutTextXY(420,1, Tss)

End

End

Until J=#27;

Serv

End;

За необхідності, можна легко модифікувати даний алгоритм, щоб він формував, наприклад, таблицю резонансних частот з будь-якою наперед заданою точністю.

Отже, маючи процедуру Specter, можна вважати задачу визначення вільних частот об'єкта вирішеною.

Зупинимося на розрахунку частотних характеристик консервативних систем. Розглянемо спочатку метод, який полягає в застосуванні формул

(4)

(5)

де - імпульсна характеристика (функція ваги) каналу об'єкта, що вивчається, і - дійсна та уявна частотні характеристики.

Імпульсну характеристику вважаємо відомою, заданою у вигляді масиву її ординат.

Цей метод досить ефективний за умови, коли імпульсна характеристика швидко затухає, тобто якщо для вона практично дорівнює нулю. Тоді верхня межа інтегрування в (4), (5) замінюється на , і застосування відповідного алгоритму чисельного інтегрування вирішує задачу.

Однак, саме для консервативних систем . Інтегрувати ж в нескінченних межах на комп'ютері немає змоги.

Отже, спробуємо змоделювати на комп'ютері реакцію каналу, що вивчається, на синусоїдальний вхідний сигнал

(6)

де не співпадає ні з однією з власних частот об'єкта. Це, як відзначалося вище, цілком реальна річ. Власні частоти також вважаємо відомими. Вимога відсутності збігу з власними частотами зрозуміла, бо частотні характеристики (принаймні, одна з характеристик (4), або (5)) на резонансній частоті буде обертатися на нескінченність. Власне, на стадії розрахунку реакції на сигнал (6) ця вимога не є категоричною, оскільки при скінченному часі спостереження вихідний сигнал хоч і буде зростати за амплітудою, але до обчислювальної катастрофи (переповнення) справа може не дійти, якщо не дуже велике. Спектральний аналіз реакції на (6), звичайно, дозволить виявити в наявність компоненти з частотою , але чисельно виділити саму цю компоненту він не дозволяє. Виділити дану складову можна, якщо скористатися методом найменших квадратів, що використовується для аналізу прихованої періодичності, вважаючи, що сигнал може бути наближено представленим (апроксимованим) виразом

. (7)

Як показник якості апроксимації береться функція

. (8)

Зрозуміло, що повинне бути мінімальним.

Умови мінімізації

(9)

Підставляючи в (9) значення Е з (8), отримуємо

,

(10)

Рівняння (10) утворюють систему з лінійних рівнянь відносно коефіцієнтів формули (7).

Коефіцієнти розширеної матриці цієї системи обчислюються за формулами-

(11)

(12)

(13)

Тут - номер рядка, - номер стовпчика (він же номер невідомого). У - ому стовпчику розміщуються праві частини рівнянь системи. Розширену матрицю формує процедура MasC.

procedure MasC;

var Z, S:integer;

Wz, Ws;real;

function Integeral (V:integer; Ws, Wz:real):real;

var i:integer;

Sum:real;

function F(t:real):real;

var t:real;

begin

t:=s*Dt;

case V of

1: F:=cos(Ws*t);

2: F:=sin(Ws*t);

3: F:=My[s];

4: F:=cos(Wz*t);

5: F:=cos(Wz*t)*cos(Ws*t);

6: F:=cos(Wz*t)*sin(Ws*t);

7: F:=My[s]*cos(Wz*t);

8: F:=sin(Wz*t);

9: F:=sin(Wz*t)*cos(Ws*t);

10: F:=sin(Wz*t)*sin(Ws*t);

11: F:=My[s]*sin(Wz*t)

end

end;

begin

Sum:=(F(0)+F(L))/2;

for i:=1 to L-1 do Sum:=Sum+F(i);

Integral:=Sum*Dt

end;

begin

С[1,1]:=D;

for S:=2 to n+2 do

С[1, S]:=Integral(1, Mw[S-2], 0);

for S:=n+3 to 2*n+3 do

СЗ[1, S]:=Integral(2, Mw[S-n-3], 0);

С[1,2*n+4]:=Integral(3,0,0);

for Z:=2 to n+2 do

begin

Wz:=Mw[Z-2];

С[Z, 1]:=Integral(4,0, Wz);

for S:=2 to n+2 do

С[Z, S]:=Integral(5, Mw[S-2], Wz);

for S:=n+3 to 2*n+3 do

С[Z, S]:=Integral(6, Mw[S-n-3], Wz);

С [Z, 2*n+4]:=Integral(7,0, Wz)

end;

for Z:=n+3 to 2*n+3 do

begin

Wz:=Mw[Z-n-3];

С [Z, 1]:=Integral(8,0, Wz);

for S:=2 to n+2 do

С [Z, S]:=Integral(9, Mw[S-2], Wz);

for S:=n+3 to 2*n+3 do

С [Z, S]:=Integral(10, Mw[S-n-3], Wz);

С [Z, 2*n+4]:=Integral(11,0, Wz)

end

end;

Процес розв’язання системи, що розглядається, здійснюється за допомогою раніше згадуваної процедури SystUr, виклик якої має вигляд -

SystUr(2*n+3, С, X).

Гармоніка з частотою в реакції на вхідний сигнал (6) буде такою:

Отже,

(14)

(15)

Повторюючи обчислення для ряду значень, визначаємо потрібні нам частотні характеристики.

Досвід експлуатації розглянутого алгоритму показує його безвідмовність, ефективність і цілком прийнятну швидкодію. Достатня точність виділення гармонічних компонент має місце за умови, коли час спостереження не менше 3 найменших і 1/3 найбільшого періодів. Кількість ординат в межах заданого часу спостереження повинна бути не меншою 20. Після округлення, умови працездатності алгоритму можуть бути записані у вигляді-

(16)

Практика експлуатації алгоритму показує, що він є стійким відносно неминучих похибок у визначенні власних частот об'єкта дослідження і похибок чисельного інтегрування початкової системи рівнянь при розрахунку _. Високочастотні складові вихідного сигналу, що відповідають вищим гармонікам власних коливань об'єкта з індексами вище , надійно відфільтровуються алгоритмом як шум. А це, в свою чергу, дозволяє обмежитися відносно невеликим значенням (для об'єктів з розподіленими параметрами інакше необхідно було б брати ). Природно, що це має місце за умови, коли .