matrične jednadžbe

Prigušivanje mehaničkih vibracija

Krešimir Burazin, Zoran Tomljanović i Ivana Vuksanović

 

Odjel za matematiku
Sveučilište Josipa Jurja Strossmayera u Osijeku
Trg Lj. Gaja 6
31000, Osijek
kburazin@mathos.hr (K. Burazin)
ztomljan@mathos.hr (Z. Tomljanović)
ivuksano@mathos.hr (I. Vuksanović)

 
Sažetak:
U radu su opisane vibracije sustava koji se sastoji od jedne mase, opruge i viskoznog prigušivača, te sustava od dvije mase, tri opruge i viskoznog prigušivača. Posebno je naglašen problem pronalaska optimalnog prigušenja slobodnih vibracija s obzirom na kriterij spektralne apscise i kriterij minimizacije ukupne energije sustava. Prikazane su i prisilne vibracije harmonijskog oscilatora, s naglaskom na pojavi rezonancije. Također su dani primjeri nekih građevina kod kojih su vibracije imale katastrofalne posljedice, te primjeri građevina kod kojih su te vibracije smanjene upotrebom raznih vrsta prigušivača.

Ključne riječi:

vibracije, harmonijski oscilator, rezonancija, prigušenje, optimalno prigušenje



1Uvod

U fizici se vibracija definira kao gibanje mase ili sustava povezanih masa oko  položaja ravnoteže. Ovdje se zapravo misli na mehaničke vibracije, odnosno na vibracije dinamičkih sustava. Međutim, sam pojam vibracije možemo shvatiti i u širem kontekstu i vjerojatno ga je, na intuitivnoj razini, najlakše usvojiti kroz razne primjere vibracionog gibanja, poput kretanja klatna sata, vibriranja žica glazbenih instrumenata, vibriranja kotača automobila, vibriranja membrane zvučnika, kucanja srca, vibriranja raznih mehaničkih uređaja (motora automobila, perilice za rublje, miksera), osciliranja serijskog strujnog kruga, vibriranja raznih građevinskih struktura (nebodera, mostova), i mnogih drugih.

Kroz ove primjere možemo uočiti da je vibracija jedan od najrasprostranjenijih oblika gibanja u prirodi. Zapravo, ukoliko probate naglas pročitati ovu rečenicu, vibriranje vaših glasnica će vam omogućiti da proizvedete potrebne zvukove, dok će vibriranje bubnjića u vašem uhu omogućiti da čujete ono što ste izgovorili.

Ljudi su se počeli zanimati za vibracije već od pojave prvih glazbenih instrumenata. Prvi koji je na znanstvenoj osnovi počeo promatrati glazbu i zvukove bio je grčki filozof i matematičar Pitagora. Kasnije su se tom vrstom gibanja bavili mnogi matematičari i fizičari, između ostalih G. Galilei, I. Newton, B. Taylor, D. Bernoulli, J. D'Alambert, L. Euler, G. R. Kirchoff, te mnogi drugi.

Već iz gore nabrojanih primjera možemo uočiti da su neke vibracije poželjne, dok neke baš i nisu. Od poželjnih vibracija, osim onih očitih koje su bitne za sam život čovjeka, poput kucanja srca i vibriranja glasnica, izdvajamo razne uređaje koji rade na principu vibracija, poput glazbenih instrumenata, kućanskih aparata (perilica rublja, mikser, shaker), satova, građevinskih strojeva (pneumatski čekić,  žaba, raketa), masažera i mnogih drugih. Nažalost, većina vibracijskih gibanja je ipak nepoželjna, i većina inženjera i znanstvenika koji se bave vibracijama zapravo pokušavaju smanjiti neželjene vibracije, odnosno bave se  prigušenjem vibracija. Neželjene vibracije se javljaju u velikom broju mehaničkih uređaja, poput turbina, pumpi, kompresora, aviona i automobila. Prevelike vibracije unutar takvih sustava proizvode neželjene zvukove, mogu dovesti do oštećenja dijelova sustava, te na kraju i do kvara, dok kod ljudi koji se nalaze u blizini izazivaju neugodu, umor i općenito gubitak efikasnosti.

Manje je poznato da i velike građevine poput mostova, nebodera i katedrala također vibriraju. Štoviše, te vibracije znaju biti prilično velike pa i uzrokovati oštećenja ili čak rušenje građevine (u zadnjem poglavlju je dano nekoliko primjera vibriranja građevina). Pri tome je izuzetno opasna pojava tzv.  rezonancije kada može doći do velikih vibracija iako sile koje djeluju na sustav nisu velike po apsolutnoj vrijednosti. Stoga se na razne načine pokušavaju smanjiti vibracije, najčešće korištenjem raznih prigušivača koji onda te vibracije svode na prihvatljivu razinu.

Matematički modeli kojima se opisuje vibracijsko gibanje obično vode na parcijalne diferencijalne jednadžbe koje su dosta složene za proučavanje, pa problemi prigušenja vibracija predstavljaju veliki izazov za znanstvenike i inženjere. Međutim, mnoge fenomene koji se pojavljuju u ovakvim sustavima možemo primijetiti već i u najjednostavnijem jednodimenzionalnom modelu harmonijskog oscilatora. U radu ćemo prikazati upravo takve vibracije promatrajući jednostavan fizikalni problem - sustav masa i opruga s jednim stupnjem slobode.

Primjeri koje dajemo u petom poglavlju pokazuju kako vibracije mogu imati razorno djelovanje na građevinske objekte. Također navodimo primjere građevina kod kojih je smanjen utjecaj vibracija te objašnjavamo na koji je način to učinjeno.



2Harmonijski oscilator

Za početak promotrimo sustav prikazan na Slici 2.1: radi se o sustavu jedne mase m\gt 0 pričvršćene elastičnom oprugom konstante elastičnosti k\gt 0 za fiksnu podlogu, s tim da je na masu istovremeno pričvršćen i prigušivač koji pruža otpor gibanju mase proporcionalan brzini gibanja, s konstantom proporcionalnosti (prigušenja) c\gt 0. Pretpostavimo da na sustav u trenutku t djeluje i vanjska sila f(t) u smjeru osi x, te da ne djeluju nikakve druge sile osim već navedenih (nema gravitacije, otpora zraka, trenja i sl.). Također pretpostavljamo da se gibanje vrši samo u smjeru osi x, te s x(t) označimo otklon tijela u trenutku t od ravnotežnog položaja. Cilj nam je pronaći položaj tijela u svakom trenutku t, odnosno pronaći funkciju pomaka x:\langle 0,\infty \rangle \rightarrow {\textbf{R}} na osnovu poznatih sila koje dijeluju na tijelo i poznatog  stanja tijela u trenutku t=0 (početni položaj i početna brzina tijela).
Img
Slika 1: Sustav mase, opruge i prigušivača.
Najprije ćemo iz Newtonovog drugog zakona: \displaystyle\ F=ma=m\frac{dv}{dt}=m\frac{d^{2}x}{dt^{2}}, gdje je a akceleracija, v brzina tijela, a F zbroj svih sila koje djeluju na masu m, odrediti jednadžbu gibanja tijela. Naime, sile koje u trenutku t djeluju na masu su
\bullet {sila opruge}, koja je prema Hookeovom zakonu proporcionalna otklonu tijela od položaja mirovanja, odnosno jednaka -kx(t);
\bullet {sila prigušenja} jednaka - cx'(t);
\bullet {vanjska sila} f(t).
Stoga drugi Newtonov zakon postaje
mx''(t)=-kx(t)-cx'(t)+f(t)\,,
odnosno dobivamo jednadžbu harmonijskog oscilatora sa prigušenjem:
(1)
mx''(t)+cx'(t)+kx(t)=f(t)\,.
To je diferencijalna jednadžba za nepoznatu funkciju x i ukoliko ju podijelimo s m\gt 0, dobivamo ekvivalentnu jednadžbu
(2)
x''(t)+ \bar{c} x'(t) +\omega_{0}^{2}x(t)=\frac{f(t)}{m},
gdje su \displaystyle\ \bar{c}=\frac{c}{m} i \displaystyle\ \omega_{0}^{2}=\frac{k}{m}, \omega_{0}\gt 0.

Taj tip diferencijalne jednadžbe naziva se linearna diferencijalna jednadžba drugog reda s konstantnim koeficijentima i dobro je poznato kako se rješava (vidi [1] za više detalja). Preciznije, ta jednadžba ima beskonačno mnogo rješenja koja se mogu eksplicitno zapisati, s tim da oblik rješenja ovisi o rješenjima pripadne karakteristične jednadžbe
\lambda^{2}+\bar{c}\lambda + \omega_{0}^{2}=0\,,
koja su dana s
(3)
\lambda_{1,2} = \frac{-\bar{c} \pm \sqrt{\bar{c}^{2}-4\omega_{0}^{2}}}{2}\,.
Razlikujemo tri moguća slučaja:
A) \lambda_{1}, \lambda_{2}\in {\textbf{R}} i \lambda_{1}\not= \lambda_{2}, tj.  \bar{c}\gt 2\omega_{0}, odnosno c\gt 2\sqrt{km}, i tada je
x(t)= x_{h}(t) + x_{p}(t)= C_{1}\text{e}^{\lambda_{1} t} + C_{2}\text{e}^{\lambda_{2} t} + x_{p}(t)\,,
gdje su C_{1}, C_{2} \in {\textbf{R}} proizvoljni, dok je x_{p} partikularno rješenje jednadžbe (2), tj.  jedno (bilo koje) rješenje te jednadžbe (ono ovisi o f i ovdje nećemo previše govoriti o metodama za pronalazak partikularnog rješenja; vidi [1]).

B) \lambda_{1},\lambda_{2} \in {\textbf{R}} i \lambda_{1}= \lambda_{2}, tj.  \bar{c}=2\omega_{0}, odnosno c=2\sqrt{km}, i tada je
x(t)= x_{h}(t) + x_{p}(t)= C_{1}\text{e}^{\lambda_{1} t} + C_{2}t\text{e}^{\lambda_{1} t} + x_{p}(t)\,,
gdje su C_{1}, C_{2} \in {\textbf{R}} proizvoljni, a x_{p} partikularno rješenje jednadžbe (2).

C) \lambda_{1}, \lambda_{2}\not\in {\textbf{R}}, tj.  \bar{c}\lt 2\omega_{0}, odnosno c\lt 2\sqrt{km}. U tom slučaju je \lambda_{1,2} = -\frac{\bar{c}}{2} \pm \bar{\omega} i, gdje je \bar{\omega} = \sqrt{\omega_{0}^{2}- \bar{c}^{2}/4}\in {\textbf{R}}\setminus \lbrace 0\rbrace, i tada je
x(t)= x_{h}(t) + x_{p}(t) = \text{e}^{-\frac{\bar{c}}{2}t} \bigl(C_{1}\sin(\bar{\omega} t) +C_{2}\cos(\bar{\omega} t)\bigr) + x_{p}(t)\,,
gdje su C_{1}, C_{2} \in {\textbf{R}} proizvoljni, a x_{p} partikularno rješenje jednadžbe (2). Lako se može provjeriti da u ovome slučaju rješenje možemo zapisati u obliku
x(t)=A\text{e}^{-\frac{\bar{c}}{2}t} \sin(\bar{\omega} t + \varphi) + x_{p}(t)\,,
gdje su A\in\textbf{R} i \varphi\in [-\pi, \pi\rangle određeni s
A=\sqrt{C_{1}^{2}+C_{2}^{2}}\,,\quad \cos \varphi = \frac{C_{1}}{\sqrt{C_{1}^{2}+C_{2}^{2}}}\,,\quad \sin \varphi = \frac{C_{2}}{\sqrt{C_{1}^{2}+C_{2}^{2}}}\,.
Napomenimo još da je za f\equiv 0 u sva tri slučaja trivijalno partikularno rješenje dano s x_{p}\equiv 0.

Uočimo da se u svakom od gore navedenih slučajeva u formuli za rješenje pojavljuju dvije proizvoljne realne konstante C_{1} i C_{2}. To ujedno znači da ta jednadžba ima beskonačno mnogo rješenja, odnosno da funkciju pomaka x nije moguće odrediti samo iz jednadžbe. Fizikalno gledano to je i očekivano, jer da bi se poznavao pomak u trenutku t nije dovoljno imati informacije samo o silama koje dijeluju na masu m. Potrebno je poznavati i stanje u kojem se masa nalazi u početnom trenutku, odnosno početni položaj x(0) i početnu brzinu x'(0) mase m.

U sustavu ne moraju biti prisutne sve navedene sile, pa tako u slučaju kada nema djelovanja vanjske sile (f\equiv 0) govorimo o slobodnim vibracijama, dok u suprotnom govorimo o prisilnim vibracijama. Isto tako, ukoliko u sustavu nema prigušivača (c =0) onda kažemo da je sustav neprigušen, a u suprotnom da je prigušen. U nastavku ćemo razmotriti slobodne vibracije, dok će o prisilnim vibracijama biti više riječi u zadnjem poglavlju.

2.1Slobodne vibracije neprigušenog sustava


Pretpostavimo najprije da na sustav ne djeluju nikakve vanjske sile i da nema prigušenja, tj. da je f\equiv 0 i c=0. Tada jednadžba (2) glasi x''(t)+ \omega_{0}^{2}x(t)=0, pripadna karakteristična jednadžba je \lambda^{2}+\omega_{0}^{2} =0, a njena rješenja su \lambda_{1,2}=\pm \omega_{0}i. Ovo je očito slučaj C, pa je rješenje dano s
(4)
x(t)=A\sin (\omega_{0}t +\varphi)\,,
gdje se konstante A i \varphi lako odrede iz početnog pomaka x(0)=x_{0} i početne brzine x'(0)=v_{0}:
(5)
A=\sqrt{x_{0}^{2}+\Big(\frac{v_{0}}{\omega_{0}}\Big)^{2}},\quad \cos\varphi=\frac{x_{0}\omega_{0}}{\sqrt{x_{0}^{2}\omega_{0}^{2}+v_{0}^{2}}}\,,\quad \sin\varphi=\frac{v_{0}}{\sqrt{x_{0}^{2}\omega_{0}^{2}+v_{0}^{2}}}\,.
Primijetimo da funkcija x poprima vrijednosti u segmentu [-A, A], te da je periodična, s temeljnim periodom \displaystyle\ T=\frac{2\pi}{\omega_{0}}. Ovakvo gibanje naziva se harmonijsko gibanje; A nazivamo amplituda, \varphi fazni pomak, a \omega_{0} vlastita frekvencija sustava. Na Slici 2.2 prikazano je jedno ovakvo gibanje za funkciju x(t)=1.5 \sin(4 t-3).
Img
Slika 2: Graf funkcije x(t)=1.5 \sin(4 t-3).


2.2Slobodne vibracije prigušenog sustava

Ukoliko je u sustavu prisutno prigušenje c\gt 0, a na njega ne djeluje nikakva vanjska sila, onda je desna strana jednadžbe gibanja (2) jednaka nuli, pa slučajevi A, B, C vode na sljedeće:

A) Ako je c\gt 2\sqrt{km} (veliko prigušenje), onda je
x(t)=C_{1}\text{e}^{\lambda_{1} t} + C_{2}\text{e}^{\lambda_{2} t}\,,
gdje su \lambda_{1,2} dani s (3), a C_{1} i C_{2} se lako izračunaju iz početnih uvjeta:
C_{1}=\frac{2v_{0}+\bar{c} x_{0} + 2x_{0}q}{4q}\,,\quad C_{2}=\frac{-2v_{0}-\bar{c} x_{0} + 2x_{0}q}{4q}\,, \quad \text{za } q=\sqrt{ \frac{\bar{c}^{2}}{4}-\omega_{0}^{2}}\,.

Uočimo da su \lambda_{1,2} \lt 0 pa je x padajuća funkcija koja teži prema nuli za velika vremena t. Ovakav oblik gibanja nazivamo pregušene oscilacije i prikazan je na Slici 2.3 zelenom bojom.

B) Slučaj c=2\sqrt{km} nazivamo kritično prigušenje i tada je (nakon što se izračunaju C_{1} i C_{2} iz početnih uvjeta)
x(t)=x_{0}\text{e}^{-\sqrt{\frac{k}{m}}\, t} + \Bigl(v_{0}+x_{0}\sqrt{\frac{k}{m}}\Bigr) t\text{e}^{-\sqrt{\frac{k}{m}}\, t}\,.
Ovakvo gibanje tijela nazivamo kritično prigušenje i jedno takvo prigušenje prikazano je na Slici 2.3 crvenom bojom.
C) Konačno, slučaj c\lt 2\sqrt{km} vodi na
x(t)=A\text{e}^{-\frac{\bar{c}}{2}t} \sin(\bar{\omega} t + \varphi)\,,
gdje je \bar{\omega} = \sqrt{\omega_{0}^{2}- \bar{c}^{2}/4}\in {\textbf{R}}\setminus \lbrace 0\rbrace, a konstante A i \varphi su određene s
A=\sqrt{x_{0}^{2}+\Big(\frac{\bar{c} x_{0}+2v_{0}}{2\bar{\omega}} \Big)^{2}}\,,\quad \cos \varphi = \frac{x_{0}}{A}\,,\quad \sin \varphi = \frac{\bar{c} x_{0}+2v_{0}}{2\bar{\omega} A}\,.
Ovakav slučaj gibanja nazivamo malo prigušenje. Veličina \displaystyle\ Ae^{-\frac{\bar{c}}{2}t} predstavlja amplitudu osciliranja i ovisi o vremenu, a \varphi fazni pomak. Uočavamo da ovo gibanje nije periodičko jer amplitude nisu konstante, već se tijekom vremena smanjuju, tj. kada vrijeme teži beskonačnosti, pomak tijela od ravnotežnog položaja teži nuli. Stoga se u ovom slučaju veličina \displaystyle\ T=\frac{2\pi}{\bar{\omega}} naziva prividnim periodom gušenja. Graf ovakvog gibanja je prikazan na Slici 2.3 plavom bojom.
Img
Slika 3: Slobodne vibracije prigušenog sustava za m=1, k=4, x_{0}=1 i v_{0}=1.


Kritično prigušenje se naziva i optimalnim prigušenjem jer odgovara prigušenju koje (u nekom smislu) optimalno prigušuje sustav. Optimalnost prigušenja bit će detaljnije diskutirana u nastavku teksta. Stoga prigušenja koja odgovaraju parametrima manjim od optimalnog nazivamo malim prigušenjima, dok veliko prigušenje odgovara parametrima koji su veći od kritičnog prigušenja.

3Sustav dvije mase

S ciljem poopćenja harmonijskog oscilatora na problem titranja više masa, promotrit ćemo sustav dviju masa koje titraju. Sustav je opisan na Slici 4 pri čemu su mase m_{1} i m_{2} međusobno povezane oprugom elastičnosti k_{2}, dok su za fiksnu podlogu povezani pripadnim oprugama elastičnosti k_{1} i k_{3}. Pretpostavit ćemo da je na prvu masu pričvršćen i prigušivač s konstantom prigušenja c, te da na mase ne djeluju dodatne sile osim navedenih (nema vanjske sile, tj. proučavamo slobodne oscilacije. Nadalje, pretpostavimo da se gibanje vrši samo u jednom smjeru te redom s x_{1}(t) i x_{2}(t) označimo otklone masa m_{1} i m_{2} od položaja ravnoteže.
Img
Slika 4: Sustav dvije mase, tri opruge i jednog prigušivača.


Da bismo opisali gibanje masa primijetimo najprije da su promjene duljina opruga k_{1}, k_{2} i k_{3} redom jednake x_{1}, x_{2}-x_{1} i -x_{2}. Sada slično kao i u slučaju harmonijskog oscilatora, uzevši u obzir drugi Newtonov zakon i Hookeov zakon, dobivamo sustav jednadžbi koji opisuje gibanje prve i druge mase:
\begin{align*} m_{1} x_{1}''(t)+c\,x_{1}'(t)+(k_{1}+k_{2})x_{1}(t)-k_{2}x_{2}(t)&=0\,,\\ m_{2} x_{2}''(t) - k_{2} x_{1}(t)+(k_{2}+k_{3})x_{2}(t)&=0\,. \end{align*}
To je sustav dvije obične diferencijalne jednadžbe s nepoznanicama x_{1} i x_{2}. Radi jednostavnije analize taj sustav ćemo zapisati u matričnom obliku:
(6)
\textbf{M}\mathsf{x}''(t)+\textbf{C}\mathsf{x}'(t)+\textbf{K}\mathsf{x}(t)=0,
gdje su nepoznanice u novom vektoru \mathsf{x}(t)=\left[ \begin{array}{c} x_{1}(t) \\ x_{2}(t) \\ \end{array} \right], a matrice su zadane s
(7)
\textbf{M}=\begin{bmatrix}m_{1}&0\\\ 0&m_{2}\end{bmatrix},\quad \textbf{K}=\begin{bmatrix}k_{1}+k_{2}&-k_{2} \\\ -k_{2}&k_{2}+k_{3} \end{bmatrix},\quad \textbf{C}=\begin{bmatrix}c&0\\\ 0&0\end{bmatrix}.
Direktnim množenjem lako se može provjeriti da je polazni sustav uistinu ekvivalentan s (6). Matrica \textbf{M} obično se naziva matrica mase, \textbf{K} matrica krutosti, dok matricu \textbf{C} nazivamo matricom prigušenja.

Sustavu (6) dodaju se i jednakosti koje određuju stanje tijela u trenutku t=0, odnosno početni položaj i početna brzina tij\'ela:
(8)
\mathsf{x}(0)=\mathsf{x}_{0} \quad \text{ i}\quad \mathsf{x}'(0)=\mathsf{v}_{0} .

3.1Slobodne vibracije neprigušenog sustava dvije mase

U početku promotrimo neprigušeni slučaj, odnosno slučaj kada je c=0. Kako bismo povezali ovaj slučaj sa slučajem harmonijskog oscilatora, moramo sustav (6) zapisati u pogodnijem obliku. U tu svrhu potrebno je uvesti pojam svojstvenih vrijednosti i svojstvenih vektora matričnog para: za vektor {\boldsymbol{\phi}} kažemo da je svojstveni vektor matričnog para (\textbf{A}, \textbf{B}) kojem pripada svojstvena vrijednost \mu \in {\textbf{C}} ako vrijedi \textbf{A}{\boldsymbol{\phi}} = \mu \textbf{B} {\boldsymbol{\phi}}, {\boldsymbol{\phi}} \neq \mathsf{0}. Poznato je da se svojstvene vrijednosti para mogu izračunati kao nultočke karakterističnog polinoma p(\mu)=\det (\textbf{A}-\mu \textbf{B}) (vidi [8]). Uočimo da ćemo za \textbf{B}=\textbf{I} dobiti standardni svojstveni problem.

Kod titrajnih problema vrlo važnu ulogu igraju svojstvene vrijednosti matričnog para (\textbf{K}, \textbf{M}), odnosno zanima nas svojstveni problem
(9)
\textbf{K}{\boldsymbol{\phi}}=\mu \textbf{M}{\boldsymbol{\phi}} .
Za matrice \textbf{M} i \textbf{K} dane u (7), lako se može provjeriti da karakteristični polinom glasi
p(\mu)= \det (\textbf{K}-\mu \textbf{M}) = m_{1} m_{2} \mu^{2} -( k_{2} m_{1} + k_{3} m_{1} + k_{1} m_{2} + k_{2} m_{2}) \mu +k_{1} k_{2} + k_{1} k_{3} + k_{2} k_{3},
pa se računanje svojstvenih vrijednosti svodi na rješavanje kvadratne jednadžbe, nakon čega se svojstveni vektori mogu odrediti rješavanjem sustava dvije jednadžbe s dvije nepoznanice, odnosno sustava (9). Lako se može vidjeti da su te svojstvene vijednosti pozitivne pa ih se može zapisati kao \mu_{i}=\omega_{i}^{2}, \omega_{i} \gt 0, i=1,2. Štoviše, pripadni svojstveni vektori {\boldsymbol{\phi}}_{1} i {\boldsymbol{\phi}}_{2} se mogu odabrati tako da vrijedi
(10)
{\boldsymbol{\Phi}}^{T} \textbf{K} {\boldsymbol{\Phi}}={\boldsymbol{\Omega}}^{2},\quad {\boldsymbol{\Phi}}^{T} \textbf{M} {\boldsymbol{\Phi}}={\textbf{I}},
gdje je
(11)
{\boldsymbol{\Phi}}=\begin{bmatrix}{\boldsymbol{\phi}}_{1}&{\boldsymbol{\phi}}_{2}\end{bmatrix}\quad \text{i} \quad {\boldsymbol{\Omega}}=\begin{bmatrix}\omega_{1}&0\\\ 0&\omega_{2}\end{bmatrix}.

Gore navedene tvrdnje vrijede i u nekim općenitijim situacijama (sustav više masa) a posljedica su toga što su \textbf{M} i \textbf{K} simetrične pozitivno definitne matrice.

Naime, za proizvoljnu simetričnu pozitivnu matricu \textbf{B} postoji matrica \textbf{C} takva da je \textbf{B}=\textbf{C}^{2}, takvu matricu \textbf{C} zovemo korijenom matrice \textbf{B} i označavamo ju s \textbf{B}^{\frac{1}{2}}. Budući da je matrica \textbf{M} simetrična pozitivno definitna tada postoji simetrična pozitivno definitna matrica \textbf{M}^{\frac{1}{2}}. Tada je matrica \textbf{K}_{1}=\textbf{M}^{-\frac{1}{2}} \textbf{K} \textbf{M}^{-\frac{1}{2}} također simetrična pozitivno definitna, stoga za nju postoji realna ortogonalna matrica \textbf{Q} i dijagonalna matrica {\boldsymbol{\Lambda}} tako da je
(12)
\textbf{Q}^{T} \textbf{K}_{1}\textbf{Q}={\boldsymbol{\Lambda}},
gdje je {\boldsymbol{\Lambda}} matrica koja na dijagonali sadrži upravo svojstvene vrijednosti matrice K_{1} koje su pozitivne zbog pozitivne definitnosti matrice \textbf{K}_{1}.

Sada za matricu {\boldsymbol{\Phi}} definiranu s {\boldsymbol{\Phi}}:=\textbf{M}^{-\frac{1}{2}}\textbf{Q} direktnim uvrštavanjem možemo provjeriti da su zadovoljene jednakosti (10). Zaista, koristeći simetričnost matrice \textbf{M}^{\frac{1}{2}} i ortogonalnost matrice \textbf{Q} slijedi
{\boldsymbol{\Phi}}^{T} \textbf{M} {\boldsymbol{\Phi}}=\textbf{Q}^{T}\textbf{M}^{-\frac{1}{2}}\textbf{M} \textbf{M}^{-\frac{1}{2}}\textbf{Q}=\textbf{Q}^{T}\textbf{Q}={\textbf{I}},
dok iz (12) dobivamo
{\boldsymbol{\Phi}}^{T} \textbf{K} {\boldsymbol{\Phi}}=\textbf{Q}^{T}\textbf{M}^{-\frac{1}{2}}\textbf{K} \textbf{M}^{-\frac{1}{2}}\textbf{Q} =\textbf{Q}^{T} \textbf{K}_{1}\textbf{Q}={\boldsymbol{\Lambda}}.
Budući da je {\boldsymbol{\Lambda}} dijagonalna matrica s pozitivnim elementima na dijagonali, to ju možemo zapisati kao {\boldsymbol{\Omega}}^{2}. Više detalja vezano uz navedene tvrdnje može se pronaći u [7, 8].

Navedene svojstvene vrijednosti i svojstveni vektori imaju veliku važnost u analizi titranja promatranog sustava. Naime, \omega_{1} i \omega_{2} se nazivaju vlastite ili svojstvene frekvencije neprigušenog sustava, a svojstveni vektori {\boldsymbol{\phi}}_{1}, {\boldsymbol{\phi}}_{2} u fizici se često zovu prvim, odnosno drugim modom titrajnog sustava, i možemo reći da oni određuju smjer gibanja pripadnih masa.

Nadalje, ako uzmemo supstituciju \mathsf{\hat{x}(t)}=\Phi^{-1} \mathsf{x}(t) u sustavu (6) (za \textbf{C}=\textbf{0}) te pomnožimo jednadžbu s lijeva s {\boldsymbol{\Phi}}^{T} i iskoristimo jednakosti dane u (10) dobit ćemo sustav
\begin{bmatrix} \hat{x}_{1}(t) \\\ \hat{x}_{2}(t)\end{bmatrix}+ \begin{bmatrix} \omega_{1}^{2}\hat{x}_{1} (t) \\\ \omega_{2}^{2}\hat{x}_{2} (t) \end{bmatrix}=\mathsf{0}.
Možemo primjetiti da smo neprigušeni problem zapisali tako da smo (u novim koordinatama) dobili titranja harmonijskih oscilatora s neprigušenim frekvencijama \omega_{1} i \omega_{2}. Stoga, postupajući analogno kao u slučaju neprigušenog harmonijskog oscilatora, dobivamo
\hat{x}_{i}(t)= A_{i} \sin (\omega_{i} t +\varphi_{i})\,,
za i=1,2, gdje se konstante A_{i} i \varphi_{i} odrede iz početnih uvjeta (8), odnosno uz spomenutu supstituciju imamo uvjete
\mathsf{\hat{x}}(0)={\boldsymbol{\Phi}}^{-1}\mathsf{x}_{0} \quad \text{ i}\quad \mathsf{\hat{x}}'(0)={\boldsymbol{\Phi}}^{-1}\mathsf{v}_{0} .
Odavde, koristeći formule (5) dobivamo
\begin{eqnarray*} &A_{i} =\sqrt{({\boldsymbol{\Phi}}^{-1}\mathsf{x}_{0})_{i}^{2}+\Big(\frac{({\boldsymbol{\Phi}}^{-1}\mathsf{v}_{0})_{i}} {\omega_{i}}\Big)^{2}},\quad \cos\varphi_{i}=\frac{({\boldsymbol{\Phi}}^{-1}\mathsf{x}_{0})_{i}\omega_{i}} {\sqrt{({\boldsymbol{\Phi}}^{-1}\mathsf{x}_{0})_{i}^{2}\omega_{i}^{2}+ ({\boldsymbol{\Phi}}^{-1}\mathsf{v}_{0})_{i}^{2}}}\,,\\ &\sin\varphi_{i} =\frac{ ({\boldsymbol{\Phi}}^{-1}\mathsf{v}_{0})_{i}}{\sqrt{({\boldsymbol{\Phi}}^{-1}\mathsf{x}_{0})_{i}^{2} \omega_{i}^{2}+{\boldsymbol{\Phi}}^{-1}\mathsf{v}_{0})_{i}^{2}}},\quad i=1,2\,. \end{eqnarray*}

Na kraju, ako iskoristimo da je \mathsf{x}(t)={\boldsymbol{\Phi}}\mathsf{\hat{x}}(t) dobit ćemo rješenje \mathsf{x}(t) u sljedećem obliku
\mathsf{x}(t)= \hat{x}_{1}(t)\boldsymbol{\phi}_{1}+ \hat{x}_{2}(t)\boldsymbol{\phi}_{2}.


3.2Slobodne vibracije prigušenog sustava dvije mase

Da bismo riješili problem s prigušenjem najprije ćemo naš sustav jednadžbi (6) zapisati u obliku sustava diferencijalnih jednadžbi prvog reda. Naime, uz već spomenutu supstituciju \mathsf{x}(t)={\boldsymbol{\Phi}}\mathsf{\hat{x}}(t) te dodatne supstitucije \mathsf{y_{1}}(t)= {\boldsymbol{\Omega}} \mathsf{\hat{x}}(t) i \mathsf{y_{2}}(t)=\mathsf{\hat{x}}'(t), sustav (6) možemo zapisati kao sustav
(13)
\mathsf{y}'(t)=\textbf{A} \mathsf{y}(t),
gdje su
(14)
\mathsf{y}(t)=\begin{bmatrix}\mathsf{y}_{1}(t)\\\ \mathsf{y}_{2}(t) \end{bmatrix},\quad \textbf{A}=\begin{bmatrix}\textbf{0}&{\boldsymbol{\Omega}}\\\ -{\boldsymbol{\Omega}}&{-}{\boldsymbol{\Phi}}^{T}\textbf{C}{\boldsymbol{\Phi}} \end{bmatrix},
a matrice {\boldsymbol{\Omega}} i {\boldsymbol{\Phi}} dane u (11). Uočimo da je za sustav dviju masa \textbf{A} matrica reda 4, dok je nepoznata funkcija \mathsf{y} vektorska funkcija koja za kodomenu ima \textbf{R}^{4}. Napomenimo i da problem harmonijskog oscilatora sa prigušenjem možemo riješiti na ovaj način, pri čemu bi \textbf{A} bila reda 2, a \mathsf{y} bi za kodomenu imala \textbf{R}^{2}.

Koristeći matrični račun rješenje dobivenog sustava možemo zapisati u obliku (vidi [1])
\mathsf{y}(t)=e^{\textbf{A}t}\mathsf{y}_{0},
gdje vektor \mathsf{y}_{0} sadrži informacije o početnim uvjetima, a matrica e^{\textbf{A}t} je definirana pomoću sume reda potencija \displaystyle e^{\textbf{A}t}=\sum_{j=0}^{\infty} \frac{\textbf{A}^{j}t^{j}}{j!} (vidi [1] i [8] za detalje). Iako efikasan izračun matrice e^{\textbf{A}t} prelazi okvire ovog rada, napomenimo da se ona može relativno jednostavno izračunati ukoliko je matrica \textbf{A} u nekoj od jednostavnijih formi, na primjer u Jordanovoj kanonskoj formi.

Kada se nađe \mathsf{y}, onda se lako izračuna \mathsf{x}:
\mathsf{x}(t)={\boldsymbol{\Phi}} \mathsf{\hat{x}}(t) = {\boldsymbol{\Phi}\boldsymbol{\Omega}}^{-1} \mathsf{y_{1}}(t)\,.



4Optimalno prigušenje

U ovome poglavlju želimo analizirati optimalno prigušenje sustava dvije mase te harmonijskog oscilatora (iako mnoge tvrdnje koje ćemo ovdje izreći vrijede i za općenitije sustave više masa). Jasno, najprije moramo precizirati što mislimo pod pojmom optimalno prigušenje. Na intuitivnoj razini možemo to ovako zamišljati: vibracije u našem sustavu su neželjene, te za dane mase i koeficijente elastičnosti želimo odrediti najbolje prigušenje c koje će osigurati najbolje (optimalno) smirivanje svih komponenti (x_{1} i x_{2}). Na primjer, na Slici 1 vidimo titranje prigušenog harmonijskog oscilatora koje odgovara raznim prigušenjima, te se možemo pitati koje prigušenje će nam biti najbolje (optimalno). Da bismo mogli odgovoriti na to pitanje najprije trebamo precizirati kriterije optimalnosti koji su u skladu s našim intuitivnim poimanjem smirivanja svih komponenti sustava, te konačno i pronaći optimalno prigušenje s obzirom na dani kriterij.

U primjeni postoji više kriterija, a mi ćemo navesti dva kriterija koja se često koriste (u obadva matrica \textbf{A} dana s (14) ima važnu ulogu):
\bullet [i)] kriterij spektralne apscise
\bullet [ii)] kriterij minimizacije ukupne energije sustava.


Kriterij spektralne apscise za optimalno prigušenje c uzima ono koje minimizira funkciju
(15)
\max_{k}{\mathrm{Re}} \lambda_{k}\,,
gdje su \lambda_{k} svojstvene vrijednosti matrice \textbf{A}. Uočimo da matrica \textbf{A} zaista ovisi o matrici \textbf{C}, odnosno o prigušenju c, pa stoga i njezine svojstvene vrijednosti možemo shvatiti kao funkcije prigušenja c. To znači da je s (15) zaista dobro definirana funkcija koja ovisi o prigušenju c, pa ima smisla tražiti njezin minimum.

Geometrijski gledano, ovaj kriterij nam kaže da je optimalno prigušenje upravo ono za koje su svojstvene vrijednosti što više lijevo u kompleksnoj ravnini.

U ovome trenutku ostavljamo nejasnom fizikalnu interpretaciju ovog uvjeta optimalnosti. Naime, iz gore navedenog nije jasno za što je takav kriterij odabran i u kojem smislu je takvo prigušenje optimalno. To ćemo pojasniti na primjeru harmonijskog oscilatora u sljedećem potpoglavlju.

Kriterij minimizacije ukupne energije sustava za optimalno prigušenje c uzima ono koje minimizira ukupnu energiju sustava definiranu sa
(16)
\int_{0}^{\infty} E(t;\mathsf{y}_{0})\, dt\,,
gdje je E(t;\mathsf{y}_{0}) ukupna energija sustava u trenutku t, odnosno zbroj potencijalne i kinetičke energije u trenutku t:
E(t; \mathsf{y}_{0})=\frac{1}{2}\dot{\mathsf{x}}(t)^{T} \textbf{M} \dot{\mathsf{x}}(t)+\frac{1}{2}\mathsf{x}(t)^{T} \textbf{K} \mathsf{x}(t)\,,
a \mathsf{y}_{0} sadrži informacije o početnim uvjetima.

Kriterij minimizacije ukupne energije intuitivno je jasan i fizikalno odgovara općenitoj ideji minimizacije energije koja se pojavljuje u raznim problemima sa fizikalnom pozadinom.

Budući da totalna energija ovisi o početnim uvjetima sadržanim u vektoru \mathsf{y}_{0}, uzimanjem usrednjenja po svim jediničnim početnim uvjetima može se pokazati (vidi npr. [5, 8]) da se ovaj kriterij svodi na minimizaciju traga matrice \textbf{X} po svim prigušenjima c, pri čemu je \textbf{X} rješenje matrične jednadžbe
(17)
\textbf{A} \textbf{X}+\textbf{X} \textbf{A}^{T}=-\textbf{I},
a matrica \textbf{A}, kao i prije, dana s (14).

Ta matrična jednadžba naziva se Ljapunovljeva jednadžba te ima veliku primjenu u teoriji upravljanja. Također uočimo da je u slučaju sustava dvije mase matrica \textbf{X} reda 4, dok je u slučaju harmonijskog oscilatora reda 2. Za matricu \textbf{A} koja ima svojstvene vrijednosti u lijevoj poluravnini kompleksne ravnine, može se pokazati da postoji jedinstveno rješenje Ljapunovljeve jednadžbe (sjetimo se da smo i kod spektralne apscise svojstvene vrijednosti pomicali što više lijevo).

Napomenimo da za ovaj kriterij postoji teorijski rezultat o optimalnom prigušenju pri čemu promatramo općenitu matricu prigušenja \textbf{C} (ne nužno oblika kao u (7)). Naime, ako bismo u jednadžbi (6) tražili optimalan \textbf{C}\in \textbf{R}^{2\times 2} za koji je \mathrm{trace}{\textbf{X}(\textbf{C})} minimalan, onda je poznato (npr. vidi [5]) da se minimum postiže za
(18)
\textbf{C}_{opt} =2 \textbf{M}^{1/2}\sqrt{\textbf{M}^{-1/2}\textbf{K}\textbf{M}^{-1/2}}\textbf{M}^{1/2},
koji se stoga često naziva kritičnim ili optimalnim prigušenjem. Istaknimo da navedeni razultat vezan uz \textbf{C}_{opt} vrijedi i općenito za sustav određen matricama n-tog reda. S druge strane, zbog fizikalne pozadine problema matrica prigušenja \textbf{C} je obično specijalne strukture (npr. kao u (7)), pa se stoga u praksi dani minimum općenito ne može postići. Međutim, za neke posebne slučajeve, poput harmonijskog oscilatora može se postići, što ćemo u nastavku i pokazati.

4.1 Optimalno prigušenje u slučaju harmonijskog oscilatora

Kada smo analizirali gibanja harmonijskog oscilatora već smo komentirali kako izgleda titranje u ovisnosti o raznim prigušenjima (vidi Sliku 1). U ovom ćemo poglavlju pokazati koliko iznosi optimalno prigušenje u slučaju harmonijskog oscilatora za kriterij spektralne apscise i za kriterij minimizacije ukupne energije sustava.

Radi potpunosti navest ćemo kako bi glasila jednadžba (13) u slučaju harmonijskog oscilatora. Naime, matrice M, C i K bile bi jednake skalarima (stoga ih ne pišemo masnim slovima), te uz oznake iz prvog poglavlja imali bismo M=m, C=c i K=k te za \Phi i \Omega iz (11) vrijedilo bi \Phi=\frac{1}{\sqrt{m}} i \Phi^{T} K \Phi = \Omega^{2}=\omega_{0}^{2}=\frac{k}{m} (odmah vrijedi i \Phi^{T} M\Phi=1). Tada bi sustav prvog reda za harmonijski oscilator glasio \mathsf{y}'(t)=\textbf{A} \mathsf{y}(t) pri čemu su funkcije \mathsf{y} i \textbf{A} dane s
(19)
\mathsf{y}(t)=\begin{bmatrix}y_{1}(t)\\\ {y}_{2}(t) \end{bmatrix},\quad \textbf{A}=\begin{bmatrix} 0 &\omega_{0}\\\ -\omega_{0} &{-}\bar{c} \end{bmatrix},
gdje je \displaystyle\ \bar{c}=\frac{c}{m}. Lako se može provjeriti da su svojstvene vrijednosti matrice \textbf{A} jednake nultočkama \lambda_{1,2} danim s (3).

Prisjetimo se da je kod kriterija spektralne apscise optimalno prigušenje ono za koje su svojstvene vrijednosti pripadne matrice \textbf{A} najviše moguće lijevo u kompleksnoj ravnini, odnosno želimo da je realni dio svojstvenih vrijednosti negativan i što veći po apsolutnoj vrijednosti.
Img
Slika 5: Svojstvene vrijednosti za k=2, m=1


Slika 5 grafički prikazuje kako se svojstvene vrijednosti, u slučaju harmonijskog oscilatora, mijenjaju u ovisnosti o prigušenju c pri čemu se c mijenja od 0.1 do 4. Kretanje svojstvenih vrijednosti matrice \textbf{A} prikazano je tako da se prigušenje mijenja od malog prigušenja 0.1 pa do velikog prigušenja koji iznosi 4. Kretanje je ilustrirano promjenom boje svojstvenih vrijednosti tako da su svojstvene vrijednosti za prigušenje 0.1 obojane sa plavom bojom, a kako parametar prigušenja raste prema broju 4 tako se boja mijenja prema crvenoj.

Možemo zaključiti da ako je c mali i ako c raste, onda su svojstvene vrijednosti kompleksne i pomiču se lijevo u kompleksnoj ravnini sve do trenutka kad svojstvena vrijednost postane dvostruka i realna te negativna. Nakon toga svojstvene vrijednosti su realne pri čemu jedna svojstvena vrijednost raste, a druga pada. Stoga, možemo zaključiti da će svojstvene vrijednosti biti najviše lijevo upravo u trenutku kad su svojstvene vrijednosti realne i dvostruke, a to je za c_{opt}=2\sqrt{km}= 2\sqrt{2}, što odgovara slučaju kritičnog prigušenja. Upravo se zato taj slučaj i naziva kritičnim (optimalnim) prigušenjem. Svojstvene vrijednosti koje odgovaraju kritičnom prigušenju označene su crvenim krugom.

Sada je vrijeme da pogledamo fizikalnu interpretaciju kriterija spektralne apscise na primjeru harmonijskog oscilatora: svojstvene vrijednosti \lambda_{1}, \lambda_{2} matrice \textbf{A} se direktno pojavljuju u formulama za funkciju pomaka danim u poglavlju 2.2. Neovisno o tome o kolikom je prigušenju riječ, pažljiviji pogled na te formule će otkriti da na amplitudu titranja najviše utječu članovi e^{{\text{R}e}\lambda_{1} t} i e^{{\text{R}e}\lambda_{2} t}. Zbog poznatog svojstva eksponencijalne funkcije da za a\le b vrijedi e^{at}\le e^{bt}, za svaki t\ge 0, to ima smisla postaviti kriterij na način da se minimiziraju {\text{R}e}\lambda_{1} i {\text{R}e}\lambda_{2}, odnosno da se minimizira veći od njih, a to je upravo ono što daje kriterij spektralne apscise. Dobru vizualizaciju gore navedenog daje i Slika 1 gdje možemo vidjeti da je u slučaju kritičnog prigušenja sustav najbrže umiren.

Pogledajmo sada koliko bi iznosilo optimalno rješenje koristeći kriterij minimizacije ukupne energije sustava. U tom slučaju potrebno je minimizirati \mathrm{trace}{\textbf{X}} pri čemu je \textbf{X} rješenje matrične jednadžbe (17) gdje je matrica \textbf{A} dana u (19).

Ljapunovljeva jednadžba s nepoznatom matricom \textbf{X} reda n može se zapisati kao sustav n^{2} linearnih jednadžbi sa n^{2} nepoznanica. U slučaju harmonijskog oscilatora nepoznata matrica je drugog reda, tj. \textbf{X}=\begin{bmatrix}{x}_{11}&{x}_{12}\\\ x_{21}& x_{22} \end{bmatrix}. Tada direktnim množenjem matrica u Ljapunovljevoj jednadžbi (za \textbf{A} kao u (19)) te izjednačavanjem komponenti s lijeve i desne strane, dobit ćemo sustav 4 jednadžbe s 4 nepoznanice:
\begin{align*} \omega_{0} x_{12} + \omega_{0} x_{21}&=-1,\\ -\omega_{0} x_{11} - \bar{c} \, x_{12} + \omega_{0} x_{22}&=0,\\ -\omega_{0} x_{11} - \bar{c} \, x_{21} + \omega_{0} x_{22}&=0, \\ -\omega_{0} x_{12} - \omega_{0} x_{21} - 2 \bar{c} \, x_{22}&=-1. \end{align*}
Budući da se radi o strukturiranom sustavu jednadžbi, lako se može izračunati da je rješenje x_{11} = 1/\bar{c } + \bar{c} /(2 \omega_{0}^{2}), x_{21}=x_{12} =-1/(2 \omega_{0}), x_{22} = 1/{\bar{c} }. Sada, koristeći \bar{c}=c/m možemo odrediti funkciju koju je potrebno minimizirati
(20)
\mathrm{trace} (\textbf{X}(c))=x_{11}+x_{22}=\frac{2m}{c}+ \frac{c}{ 2 \omega_{0}^{2} m},\qquad \text{gdje je }c\gt 0\text{}.
Budući da se radi o glatkoj funkciji jedne varijable, standardnim postupkom minimizacije iz
\frac{d \mathrm{trace} (\textbf{X}(c))}{dc}= \frac{-2m}{c^{2}}+ \frac{1}{2 \omega_{0}^{2} m}
dobivamo stacionarnu točku c_{opt}=2\sqrt{\omega_{0}^{2} m^{2}}. Koristeći \omega_{0}^{2}=k/m slijedi c_{opt}=2\sqrt{km}, što se, očekivano, podudara s kritičnim prigušenjem iz poglavlja 2.1. Uočimo i da je funkcija c \mapsto \mathrm{trace} (\textbf{X}(c)) strogo padajuća (derivacija joj je negativna) na \langle 0, c_{opt}\rangle te strogo rastuća (derivacija joj je pozitivna) na \langle c_{opt}, +\infty\rangle, pa možemo zaključiti da je dobiveni c_{opt} točka globalnog minimuma.

Primjetimo da se dobivena optimalna prigušenja za kriterij minimizacije ukupne energije sustava i kriterij spektralne apscise podudaraju: c_{opt}=2\sqrt{km}. Nadalje, ako bismo izračunali optimalni minimum po svim prigušenjima koji je dan formulom (18) dobili bismo također c_{opt}=2\sqrt{km}. Međutim, važno je istaknuti da to općenito nije slučaj. Naime, već za slučaj dvije mase globalni minimum dan s (18) se ne postiže, odnosno pomoću prigušivača u našem modelu ne možemo postići kritično prigušenje. U sljedećem poglavlju ćemo vidjeti na primjeru sustava dvije mase i da se optimalno prigušenje dobiveno minimizacijom ukupne energije sustava općenito ne podudara s optimalnim prigušenjem koji daje kriterij spektralne apscise. Isti se fenomen može primijetiti i kod općenitih titrajnih sustava.

4.2 Optimalno prigušenje u slučaju dvije mase


Problem minimizacije funkcije t(\textbf{C}):=\mathrm{trace} \textbf{X(C)} općenito nije jednostavan minimizacijski problem, te često vodi na nekonveksni optimizacijski problem. Općenito, prilikom optimizacijskog postupka potrebno je riješiti velik broj Ljapunovljevih jednadžbi sa matricama velikih dimenzija što je vremenski i numerički vrlo zahtjevno. Rješavanje Ljapunovljeve jednadžbe jedno je od važnih tema u području numeričke linearne algebre gdje se posebna pozornost pridaje što točnijem efikasnom rješavanju te jednadžbe (za neke rezultate vidi [8, 2, 3]). Zbog spomenute složenosti problema, nećemo ulaziti u detaljne dokaze, ali ćemo ilustrirati dio osnovnih rezultata.

Uočimo da bi u slučaju titranja dvije mase matrice koje se pojavljuju u Ljapunovljevoj jednadžbi bile reda 4, što bi značilo da trebamo riješiti 16 jednadžbi sa 16 nepoznanica kako bismo odredili \mathrm{trace} \textbf{X}, pa stoga sustav jednadžbi nećemo precizno zapisivati. Međutim, istaknimo da se, u slučaju dvije mase s jednim prigušivačem c, može pokazati da {\mathrm{trace}}\textbf{X}(c) ima oblik (vidi npr. [9]) {\mathrm{trace}}\textbf{X}(c)=\frac{A}{c}+ B\,c, gdje su A i B konstante dane u terminima početnih parametara našeg sustava. Uočite da se izgled funkcije podudara sa funkcijom dobivenom u (??). Stoga, potpuno analogno kao u prošlom poglavlju može se pokazati da je optimalno prigušenje c_{opt} dano je formulom c_{opt}= \sqrt{A/B}.

Što se tiče kriterija spektralne apscise, na primjeru ćemo skicirati koliko iznosi optimalno prigušenje te ćemo ga usporediti s optimalnim prigušenjem u slučaju kriterija minimizacije ukupne energije sustava. U tu svhu promotrimo problem titranja dviju masa koji je dan na Slici 4 za sljedeće parametre sustava m_{1}=20, \, m_{2}=40, k_{1}=k_{2}=k_{3}=1, dok vrijednost parametra c varira od 0.1 do 20.

Img
Slika 6: Svojstvene vrijednosti za sustav dvije mase


Slika 6 prikazuje kretanje svojstvenih vrijednosti matrice \textbf{A} pri čemu se parametar prigušenja mijenja od malog prigušenja 0.1 pa do velikog prigušenja koji iznosi 20. Analogno kao na Slici 5, kretanje je ilustrirano promjenom boje svojstvenih vrijednosti tako da su svojstvene vrijednosti za prigušenje 0.1 obojane sa plavom bojom, a kako parametar prigušenja raste prema broju 20 tako se boja mijenja prema crvenoj. Sa slike možemo primijetiti da će optimalan c prema kriteriju spektralne apscise biti upravo kada su svojstvene vrijednosti najviše lijevo, što odgovara onom prametru c kada smo dobili svojstvene vrijednosti na slici označene crvenim krugom. Kritično prigušenje postiglo bi se u slučaju kada su sve svojstvene vrijednosti realne, a već sa slike je očito da se to ne može dogoditi.

Pokažimo na kraju da optimalna prigušenja za kriterij spektralne apscise i kriterij minimizacije ukupne energije sustava nisu jednaka. Prisjetimo se da optimalna prigušenja odgovaraju minimumima funkcija r(c):= \max_{k} \mathrm{Re} \lambda_{k} za kriterij spektralne apscise, odnosno t(c):=\mathrm{trace}{\textbf{X}(c)}, za kriterij minimizacije ukupne energije sustava. Grafovi tih funkcija (za sustav dvije mase s parametrima m_{1}=20, \, m_{2}=40, k_{1}=k_{2}=k_{3}=1, kao i prije) prikazani su na slici 7.
Img
Slika 7: Usporedba optimalnih prigušenja za različite kriterije


Iz Slike 7 se jasno vidi da im minimumi nisu isti, odnosno optimalna prigušenja za ova dva kriterija nisu ista. Preciznije, optimalno prigušenje za kriterij spektralne apscise iznosi 7.125, dok za kriterij minimizacije ukupne energije iznosi 7.442. Prisjetimo se da su kod harmonijskog oscilatora optimalna prigušenja bila jednaka, ali kao što je ovdje ilustrirano, to općenito nije slučaj.

U ovom poglavlju su ilustrirane sličnosti, ali i bitne razlike između problema optimalnog prigušenja kod titranja harmonijskog oscilatora i problema titranja dvije mase. Slični fenomeni se pojavljuju i kod titranja složenijih mehaničkih sustava, koji se često mogu diskretizirati i zapisati sustavom (6), ali sa matricama sustava puno veće dimenzije. Tada složenost problema znatno raste, a za određivanje optimalnog prigušenja korisno je i razumijevanje jednostavnijih problema poput onih prikazanih u ovom radu.

U kontekstu optimizacije prigušivača napomenimo da se možemo pitati gdje postaviti prigušivač tako da je sustav optimalno prigušen (u literaturi to zovemo optimizacija položaja prigušivača). U našem slučaju dvije mase, to znači treba li prigušivač postaviti na prvu ili drugu masu. Općenito optimizacija prigušivača (optimizacija položaja i njihovih optimalnih viskoznosti) rezultira nekonveksnim optimizacijskim problemom pri čemu se efikasni optimizacijski algoritam još uvijek traži. Zbog kompleksnosti optimizacijskog problema, obično se pri određivanju optimalnog prigušenja koriste aproksimacijske tehnike (npr. vidi [3, 2]).

Neki složeniji poznati primjeri vibracija građevinskih objekata navedeni su u sljedećem poglavlju.



5Rezonancija

Do sada smo promatrali optimalno prigušenje slobodnih oscilacija harmonijskog oscilatora i sustava dvije mase. U ovome poglavlju ćemo po prvi puta promatrati sustav na koji dodatno djeluje neka vanjska sila, odnosno promatrat ćemo prisilne oscilacije harmonijskog oscilatora, s posebnim naglaskom na pojavi rezonancije. Složeniji problem prisilnih oscilacija sustava dvije mase ovdje nećemo promatrati.

Preciznije, pretpostavimo da na harmonijski oscilator iz drugog poglavlja djeluje periodična sila oblika f(t)=F\cos (\omega t), za neke F,\omega \in {\textbf{R}}^{+}. Da bi pronašli funkciju pomaka x sada je potrebno još pronaći i partikularno rješenje x_{p} jednadžbe (2). Za ovakav oblik funkcije smetnje f to se može napraviti tzv. metodom neodređenih koeficijenata, koja se sastoji u tome da partikularno rješenje pokušamo pronaći u obliku
(21)
x_{p}(t)=A_{1}\cos (\omega t)+A_{2} \sin (\omega t)\,,
pri čemu su A_{1} i A_{2} nepoznate konstante. Uvrštavanjem ovakvog x_{p} i njegovih derivacija u (2) lako se izračunaju A_{1} i A_{2}, te se tako dobije formula za rješenjenje.

Pogledajmo najprije što se događa u sustavu bez prigušenja (c=0): ukoliko je \omega\not=\omega_{0} dobije se \displaystyle\ A_{1}=\frac{F}{m(\omega_{0}^{2}- \omega^{2})}, A_{2}=0, pa je opće rješenje jednadžbe (2) jednako
x(t) = A \sin (\omega_{0} t + \varphi)+\frac{F}{m(\omega_{0}^{2}- \omega^{2})}\cos (\omega t)\,,
pri čemu A i \varphi određujemo iz početnih uvjeta, kao i prije. Uočimo da je ovdje gibanje opisano kao zbroj (superpozicija) dvaju harmonijskih gibanja s periodima \displaystyle\ T_{1}=\frac{2\pi}{\omega_{0}} i \displaystyle\ T_{2}=\frac{2\pi}{\omega}, te pripadnim amplitudama A i \displaystyle\ \frac{F}{m(\omega_{0}^{2}- \omega^{2})}. Primjetimo i da amplituda \displaystyle\ \frac{F}{(\omega_{0}^{2}- \omega^{2})} ovisi o vanjskoj sili te da može biti velika u dva slučaja: ukoliko je amplituda vanjske sile F velika, te ukoliko je \omega približno jednak \omega_{0}. U graničnom slučaju, kada je \omega_{0}=\omega, tj. kada je vlastita frekvencija sustava jednaka frekvenciji vanjske sile, dolazi do pojave poznate pod nazivom rezonancija. Tada partikularno rješenje ne možemo pronaći u obliku (21), nego ga tražimo kao
x_{p}(t)=t(A_{1} \cos (\omega t)+ A_{2} \sin (\omega t))\,.
Analogno kao i prije se izračuna A_{1}=0 i A_{2}=\frac{F}{2m\omega}, pa je pomak dan s
x(t)=A \sin (\omega_{0} t + \varphi)+\frac{Ft}{2m\omega}\sin (\omega t)\,,
pri čemu se A i \varphi izračunaju iz početnih uvjeta. I ovdje je pomak superpozicija dvaju gibanja, od kojeg je prvo harmonijsko, dok se kod drugog amplituda linearno povećava s vremenom. Stoga će za dovoljno veliki t amplituda postati prevelika, te će doći do rušenja sustava. Slika 8 prikazuje primjer prisilne vibracije pri pojavi rezonancije.
Img
Slika 8: Pojava rezonancije (m=1, k=4, x_{0}=1, v_{0}=1, F=1, \omega =1)


Rezonancija je jedan od uzroka oštećivanja raznih mehaničkih sustava, a neki poznatiji primjeri rušenja građevinskih konstrukcija kao posljedica rezonancije dani su na kraju ovog poglavlja.

Kako bismo spriječili neželjene posljedice rezonancije, u sustav na koji djeluje vanjska sila f možemo dodati prigušivače. Promotrimo sada naš sustav s vanjskom silom f(t)=F\cos (\omega t), s tim da sada u sustavu imamo prigušenje c\gt 0. Tada je on opisan nehomogenom diferencijalnom jednadžbom
(22)
x''(t)+\bar{c} x'(t)+\omega_{0}^{2}x(t)=\frac{F}{m} \cos (\omega t)\,.
Vidjeli smo da izgled rješenja ovisi o nultočkama pripadne karakteristične jednadžbe (slučajevi A,B,C iz poglavlja 2.1). U svakom od tih slučajeva partikularno rješenje možemo pronaći, kao i kod neprigušenih prisilnih vibracija, u obliku
x_{p}(t)=A_{1}\cos (\omega t)+A_{2} \sin (\omega t)\,.
I ovaj put nepoznate konstante A_{1} i A_{2} dobivamo deriviranjem i uvrštavanjem u (22), s tim da ovaj put dobivamo nešto kompliciranije izraze:
\begin{eqnarray*} &&A_{1}=\frac{F(\omega_{0}^{2}-\omega^{2})/m}{(\omega_{0}^{2}-\omega^{2})^{2}+(\bar{c} \omega)^{2}}\,, \\ &&A_{2}=\frac{F\bar{c} \omega/m}{(\omega_{0}^{2}-\omega^{2})^{2}+(\bar{c} \omega)^{2}}. \end{eqnarray*}
Prema tome, partikularno rješenje je oblika
x_{p}(t)=A_{1}\cos (\omega t)+A_{2} \sin (\omega t)=\widetilde{A} \sin (\omega t+ \psi),
gdje je
\begin{eqnarray*} \widetilde{A}=\sqrt{A_{1}^{2}+A_{2}^{2}}&=& \frac{F}{m\sqrt{(\omega_{0}^{2}-\omega^{2})^{2}+(\bar{c} \omega)^{2}}}\,,\\ \sin\psi &= &\frac{\omega_{0}^{2} - \omega^{2}}{\sqrt{(\omega_{0}^{2}-\omega^{2}) + (\bar{c}\omega)^{2}}}\,,\\ \cos\psi &= &\frac{\bar{c}\omega}{\sqrt{(\omega_{0}^{2}-\omega^{2}) + (\bar{c}\omega)^{2}}}\,. \end{eqnarray*}
Sada je funkcija pomaka (opće rješenje jednadžbe (22)) superpozicija prigušenih slobodnih vibracija x_{h} i vibracija x_{p} koje uzrokuje vanjska sila:
x=x_{h}+x_{p}= x_{h} +\widetilde{A} \sin (\omega t + \psi)\,.
Pri tome x_{h} ovisi o veličini prigušenja c, pa tako možemo imati malo, veliko ili kritično prigušenje. U svakome od ta tri slučaja smo vidjeli da pomak x_{h}(t) teži nuli kada t \to \infty i on predstavlja tzv. prijelazno rješenje. Drugi dio, x_{p}, je periodična funkcija s periodom \displaystyle\ T=\frac{2 \pi}{\omega} i ne ovisi o početnim uvjetima, već samo o vanjskoj sili. Primjer takvog ponašanja sustava dan je na Slici 9.
Img
Slika 9: Prisilne vibracije prigušenog sustava (m=1, k=4, x_{0}=1, v_{0}=1, F=1, \omega =1)


Uočimo da ukoliko je \bar{c} dovoljno malen i \omega \approx \omega_{0}, amplituda \widetilde{A} može biti vrlo velika, pa i u ovom slučaju može doći do rezonancije. Stoga se prirodno postavlja pitanje optimalnog prigušenja za prisilne oscilacije. Međutim, to je još složeniji problem od prigušenja slobodnih oscilacija i prelazi okvire ovoga rada.

Na kraju, da bismo još jednom istaknuli važnost problema neželjenih vibracija, navodimo nekoliko primjera vibriranja građevina.

Tacoma Bridge

U srpnju 1940. godine bio je završen i otvoren most u Washingtonu poznat pod nazivom Tacoma Narrows Bridge. Bio je to jedan od najvećih visećih mostova tog vremena. Osiguravala su ga užad koja su jednim krajem bila pričvršćena za tlo, a drugim za dijelove mosta koji su se najviše micali. Od prvoga dana most je počeo vertikalno oscilirati. Nakon samo tri mjeseca užad namijenjena stabiliziranju mosta su popucala tijekom oluje. Iako su zamijenjena, most je sve više oscilirao prilikom jačih vjetrova. U studenome iste godine, uz brzinu vjetra od oko 65 km/h,u jednom trenutku jedan kraj ceste bio je 8 metara viši od drugog, a u drugom pak 8 metara niži od drugog. Ubrzo nakon toga most se raspao. Uzrok ove katastrofe jest upravo rezonancija: uslijed izjednačavanja frekvencije mosta i frekvencije vjetra amplituda titranja mosta je postala prevelika, te se on srušio.

Kasnije je most je obnovljen, no sada ima nosače koji su 10 metara duboko u tlu, a ne 2.4 metra kao kod prvotnog mosta i ima četiri, a ne dvije trake. Masa novog mosta je duplo veća od starog, i još uvijek je u uporabi ([4, 10]).

Millenium bridge

Milenijski most, poznat kao i Londonski milenijski pješački most je čelični viseći pješački most preko prijeke Temze u Londonu. Otvoren je u lipnju 2000. godine, a dobio je nadimak Wobbly Bridge (Nesiguran most) nakon što su pješaci osjetili neočekivano ljuljanje u prva dva dana nakon otvorenja. Most je tada zatvoren, te se nije otvarao iduće dvije godine, sve dok ljuljanja nisu posva uklonjena. Što je bio uzrok gibanja mosta? Na dan otvaranja most je prešlo 90 000 ljudi, čija je kretnja uzrokovala male bočne oscilacije mosta. Ne samo Milenijski most, nego i bilo koji most bočne frekvencije manje od 1.3 Hz i dovoljno male mase može svjedočiti istom problemu. Što je veći broj ljudi na mostu, to je i veća amplituda vibracija. Nakon opsežne analize inženjeri su uklonili problem Milenijskog mosta ugradnjom 37 viskoznih prigušivača koji kontroliraju horizontalno gibanje, te 52 masovna prigušivača koji kontroliraju vertikalno gibanje mosta. Nakon provedenih testiranja, most je otvoren, te otada nije bilo nikakvih značajnih vibracija mosta ([11]).

Viseći most u Broughtonu

Pojava rezonancije je odgovorna i za urušavanje visećeg mosta u Broughtonu pored Manchestera, 1831. godine. Ovo se dogodilo kada je kolona vojnika marširala preko mosta, stvarajući pri tome periodičnu silu dosta velike amplitude. Frekvencija te sile bila je jednaka prirodnoj frekvenciji mosta. Tako su bile pobuđene vrlo velike vibracije i most se raspao. Upravo zbog ovoga razloga dana je naredba vojnicima da se marš prekida prilikom prelaska preko mosta ([4]).

Most Franje Tuđmana

Most Franje Tuđmana u Dubrovniku sagrađen je 2001. godine. Asimetrična konstrukcija ima 143 metra visok toranj na koji je sa svake strane pričvršćeno je 19 čeličnih konopa od kojih najveći ima duljinu 220 metara. U ožujku 2005. i 2006. godine snježne oluje su uzrokovale ekstremne vibracije tih konopa, od kojih je najduži vibrirao amplitudom do 2 metra. Stoga su stručnjaci odlučili postaviti vibracijske prigušivače. Postavljeni su vertikalno između ceste i konopa, te pričvršćeni za konop otprilike 3.5 metra iznad ceste. Rezultati su pokazali da su vibracije smanjene za faktor 10. Iskustvo pokazuje da bi ovi prigušivači trebali osigurati da se u budućnosti, čak i pod uvjetima snažnog vjetra, konoplje mosta zanjiše s malom amplitudom koja nije presudna građevinskoj sigurnosti. Izračuni pokazuju da amplituda vibriranja konoplja u najgorem slučaju neće prelaziti 15-20 cm ([12]).

Burj al Dubai i Taipei

Burj al Dubai, poznat i kao Burj Khalifa, najviša je svjetska građevina, visine 828 metara, sa 206 katova. Posebno je dizajniran kako bi se smanjio utjecaj snage vjetra na toranj. Naime, temelji Bujr al Dubai-a postavljeni su u obliku slova "Y", a zgrada se diže u nebo u nekoliko odvojenih stupova koji se uzdižu oko centralnog tornja. Taj dizajn odmiče vjetar od građ evine i sprječava ga u stvaranju vrtloga zračnih struja koje bi uzrokovale ljuljanje zgrade. Usprkos ovome strateškom dizajnu, Burj Khalifa se njiše pri vrhu amplitudom 2 metra. Poput mnogih nebodera, Bjur Khalifa za smanjivanje utjecaja vjetra koristi i tzv. strateške masivne prigušivače. To su ogromna njihala čiji oblik i veličina ovise o masi i visini svakog odvojenog tornja. Kako vjetar puše i gura zgradu u jednom smjeru, tako prigušivač klizi u suprotnom smjeru, smanjujući tako ljuljanje zgrade.
Ovakav prigušivač ima i druga po redu najviša svjetska građevina, neboder Taipei 101 u Tajvanu. Unutar nebodera, između 88. i 92. kata smješteno je ogromno njihalo (prigušivač) koje vodi tihu bitku s velikim olujama i tajfunima. Metalna kugla, težine 730 tona, blago se njiše naprijed-nazad, te tako štiti neboder od sile vjetra i osigurava udobnost njegovim stanarima ([13, 14]).

Katedrala sv. Petra u Beauvais-u

Katedrala sv. Petra u Beauvais-u u Francuskoj visoka je 153 metra i jedno je od najdragocjenijih postignuća gotičke arhitekture. Sa željom da se sagradi najviša katedrala u 13. stoljeću, graditelji su prešli granice tadašnje tehnologije. Iako je katedrala bila tada najviša na svijetu, imala je relativno tanke potporne stupove kako bi što više svjetlosti moglo ući u nju. 12 godina nakon što je sagrađena, kor katedrale se urušio, zajedno s nekoliko gornjih potpornja. Za ovu katastrofu odgovorna je rezonancija koja je nastala pod utjecajem snažnih vjetrova koji su došli s engleskog kanala. Kor je obnovljen, a dodani su i masivni potpornji kako bi stabilizirali sjevernu stranu katedrale. Iako je broj potpornja udvostručen, i danas katedrali prijete olujni udari vjetrova koji uzrokuju njihovo vibriranje i slabljenje krovnih greda. Znanstvenici diljem svijeta proučavaju katedralu i pokušavaju naći rješenje kako bi spriječili novu katastrofu ([15, 16]).




6Literatura
References
[1] M. Alić, Obične diferencijalne jednadžbe, Matematički odjel Prirodoslovno-matematičkog fakulteta, Sveučilište u Zagrebu, Zagreb, 1994.
[2] P. Benner, Z. Tomljanović, and N. Truhar, Dimension reduction for damping optimization in linear vibrating system, Journal of Applied Mathematics and Mechanics (ZAMM), 2011, 91 (3), 179-191; DOI: 10.1002/zamm.201000077.
[3] P. Benner, Z. Tomljanović, and N. Truhar, Optimal damping of selected eigenfrequencies using dimension reduction, Numerical Linear Algebra with Applications, 20 (2013), no. 1, 1–17, DOI: 10.1002/nla.833.
[4] M. Braun, Differential equations and their applications, Springer, New York, 1992.
[5] S. Cox, I. Nakić, A. Rittmann, and K. Veselić, {Lyapunov} optimization of a damped system, Systems & Control Letters 53 (2004), 187–194.
[6] J. M. Krodkiewski, Mehanical Vibrations, The University of Melbourne, 2008.
[7] N. Truhar, Numerička linearna algebra, Odjel za matematiku, Svučilišta J. J. Strosmayera u Osijeku, 2010.
[8] K. Veselić, Damped oscillations of linear systems – a mathematical introduction, Springer, 2011.
[9] K. Veselić On linear vibrational systems with one-dimensional damping, Applied Analysis 29 (1988), 1–18.
[10] http://matdl.org/failurecases/Bridge_Collapse_Cases/Tacoma_Narrows (13.7.2011.)
[11] http://en.wikipedia.org/wiki/Millennium_Bridge_(London) (9.9.2011.)
[12] http://www.empa.ch/plugin/template/empa/*/51076/ (pristup 26.7.2006.)
[13] http://www.enggpedia.com/civil-engineering-encyclopedia/megastructures/b... (15.9.2011.)
[14] http://www.studentpulse.com/articles/124/4/confusing-the-wind-the-burj-k... (2011)
[15] http://www.learn.columbia.edu/ma/htm/ms/ma_ms_bc_discuss_collapse.htm (2000.)
[16] http://archive.cyark.org/cathedral-of-beauvais-info (1.9.2006.)