mehanika

Metoda konačnih elemenata; teorija i praktična implementacija

Andrej Novak
Sveučilište u Dubrovniku
1Uvod

Metoda konačnih elemenata je jedna od najpopularnijih numeričkih metoda za rješavanje parcijalnih diferencijalnih jednadžbi. Cilj ovog rada je na jednostavan način, kroz primjere, prezentirati metodu konačnih elemenata. U prvom dijelu rada je opisana Galerkinova metoda te je kroz jednostavan primjer skrenuta pozornost na opasnosti prilikom numeričke implementacije. U središnjem dijelu rada je opisana metoda konačnih elemenata kao poseban slučaj Galerkinove metode, a kroz primjer je dano numerički prihvatljivije rješenje primjera iz prvog dijela rada, samo što je jedan rub ostavljen slobodnim kako bi se prikazala implementacija raznih rubnih uvjeta. U završnom dijelu rada su dane neke ideje o implementaciji metode konačnih elementa na problemima koji uključuju više derivacije te par komentara vrijednih za praktičnu implementaciju.

2Galerkinova metoda

U ovom odjeljku ćemo postaviti scenu za metodu konačnih elemenata. Započinjemo kratkim opisom Galerkinove metode. Neka je V Hilbertov prostor, a(\cdot , \cdot) : V \times V \to \mathbb{R} bilinearna forma koja je ograničena,
 

|a(u,v)| \leq M\Vert u\Vert _{V}\Vert v\Vert _{V} \forall u, v \in V

i V-eliptična,
 

a(v,v) \geq c \Vert v\Vert _{V}^{2} \forall v \in V

te l ograničeni linearni funkcional. Promatramo problem

(1)
u \in V, \hspace{5mm} a(u,v) = l(v) \hspace{5mm} \forall v \in V.

Prema Lax-Milgramovoj lemi, varijacijski problem (1) ima jedinstveno rješenje. Općenito govoreći, nemoguće je naći egzaktno rješenje od (1) zato što je V beskonačno dimenzionalni prostor. Prirodno je probati konstruirati aproksimaciju rješenja u konačnodimenzionalnom prostoru V_{N} \subset V. Promotrimo projekciju problema (1) na V_{N},

(2)
u_{N} \in V_{N}, \hspace{5mm} a(u_{N},v) = l(v) \hspace{5mm} \forall v \in V_{N}.

Uz pretpostavke da je bilinearna forma a ograničena, V - eliptična i l ograničeni linearni funkcional, opet primjenom Lax-Milgramove leme slijedi da (2) ima jedinstveno rješenje u_{N}.


Neka je \lbrace \phi_{i}\rbrace _{i=1}^{N} baza konačnodimenzionalnog prostora V_{N}, tada možemo zapisati

u_{N} = \sum_{j = 1}^{N} \xi_{j}\phi_{j}

za neke, trenutno nepoznate, skalare \xi_{i}, i = 1,2,..,N.
Uzmimo da v \in V_{N} budu baš \phi_{i} te uvrstimo u (2). Dobivamo

(3)
\sum_{j = 1}^{N} a( \phi_{j},\phi_{i})\xi_{j} = l(\phi_{i}) \hspace{5mm} i = 1,...,N.

Očito, prethodnu tvdnju možemo zapisati kao linearni sustav

\mathbf{K\xi = b}

gdje je vektor nepoznanica \mathbf{\xi} = (\xi_{j}) \in \mathbb{R}^{N}, \mathbf{K} = (a(\phi_{j},\phi_{i})) \in \mathbb{R}^{N\times N} matrica sustava te \mathbf{b} = (l(\xi_{i}))\in \mathbb{R}^{N} desna strana.


Općenito je aproksimacija u_{N} različita od rješenja u, stoga je prirodno tražiti približno rješenje u_{N} u potprostoru V_{N} što veće dimenzije. Za niz potprostora V_{N_{i}} takvih da vrijedi V_{N_{1}} \subset V_{N_{2}} \subset ... V_{N}, računamo pripadna rješenja u_{N_{i}}, i = 1,2,.... Ovaj postupak nazivamo Galerkinovom metodom. Prvi korak u dokazu konvergencije Galerkinove metode je važan i elegantan teorijski rezultat.

Teorem 1. (Céina lema) Neka je V Hilbertov prostor, V_{N} \subset V,a( \cdot, \cdot) ograničena, V-eliptična bilinerna forma na V te l ograničeni linearni funkcional. Neka je u \in V rješenje od (1) i u_{N}\in V_{N} Galerkinova aproksimacija (2). Tada postoji konstanta C takva da
(4)
\Vert u - u_{N}\Vert _{V} \leq C\inf_{v \in V_{N}} \Vert u - v\Vert _{V}.

Proof. Prema pretpostavci a( \cdot, \cdot) je V-eliptična, ograničena bilinearna forma pa za svaki v\in V_{N} imamo,
\begin{align*} c\Vert u - u_{N}\Vert _{V}^{2} &\leq a(u - u_{N}, u - u_{N})\\\& = a(u- u_{N}, u - v)\\\ &\leq M\Vert u - u_{N}\Vert _{V}\Vert u-v\Vert _{V}. \end{align*}
Dijeljenjem s c\Vert u - u_{N}\Vert slijedi tvrdnja,
\Vert u - u_{N}\Vert _{V}^{2} \leq C\Vert u - u_{N}\Vert _{V}\Vert u-v\Vert _{V}
gdje je C = M/c.
\ \blacksquare

U slučaju kada je a(\cdot,\cdot) simetrična bilinearna forma ona definira skalarni produkt na V i njime induciranu normu \Vert v\Vert _{a} = \sqrt{a(v,v)}. Ako sada oduzimemo (2) od (1) zbog linearnosti slijedi

a(u-u_{N},v) = 0 \hspace{5mm} \forall v\in V_{N}.

Vidimo da je s obzirom na novi skalarni produkt pogreška u - u_{N} ortogonalna na potprostor V_{N}. Drugim riječima, u_{N} je ortogonalna projekcija egzaktnog rješenja u na potprostor V_{N}.

Korolar 2. Uz pretpostavke prethodnog teorema pretpostavimo još da je dan niz potprostora V_{N_{1}} \subset V_{N_{2}} \subset ... takvih da je \overline{\bigcup_{i}V_{N_{i}}} = V. Tada Galerkinova metoda konvergira,
\Vert u - u_{N_{i}}\Vert \longrightarrow 0 \hspace{4mm} \text{ kad }i \longrightarrow \infty.

Proof. Zbog pretpostavke gustoće postoji niz v_{i} \in V_{N_{i}}, i \geq 1, takav da
\Vert u - v_{i}\Vert _{V}\longrightarrow 0 \hspace{4mm} \text{ kad } i \longrightarrow \infty.
Primjenom Céine leme slijedi dokaz tvrdnje,
\Vert u - u_{N_{i}}\Vert _{V} \leq C\Vert u-v_{i}\Vert _{V}.
\ \blacksquare

Promotrimo primjenu Galerkinove metode na jednostavnom primjeru.

Primjer 3. Zadan je rubni problem
(5)
\begin{cases}-u'' + u = f \text{ na } \langle0,1\rangle\\ u(0) = u(1) = 0. \end{cases}
Množenjem s test funkcijom v \in V = H_{0}^{1}(\langle0,1\rangle), integriranjem po cijeloj domeni te primjenom parcijalne integracije slijedi slaba formulacija
u\in V, \hspace{5mm} \int_{0}^{1}(u'v' + uv)dx = \int_{0}^{1}fvdx, \hspace{5mm} \forall v\in V.
Ako definiramo a(u,v) = \int_{0}^{1}(u'v' + uv)dx te l(v) = \int_{0}^{1}fvdx iz Lax-Milgramove leme slijedi jedinstvenost rješenja ovako zapisane zadaće.
Moramo odabrati konačnodimenzionalni potprostor od V u kojem želimo tražiti aproksimaciju rješenja. S obzirom na to da je rješenje u rubovima nula jedan jednostavan potprostor je
V_{N} = span\lbrace x^{i}(1-x), i = 1,2,...,N\rbrace .
Onda je
u_{N} = \sum_{j = 1}^{N} \xi_{j}x^{j}(1-x).
Nepoznate koeficijente \lbrace \xi_{j}\rbrace _{j = 1}^{N} odredit ćemo iz Galerkinovih jednadžbi
\int_{0}^{1}(u_{N}'v' + u_{N}v)dx = \int_{0}^{1} fv dx\hspace{5mm} \forall v\in V_{N},
uzimanjem funkcija v = x^{i}(1-x),i= 1,2,...,N za test funkcije. Dolazimo do sustava linearnih jednadžbi
\mathbf{K\xi = b},
gdje je {\mathbf{\xi} }= (\xi_{1},...,\xi_{N})^{\tau} vektor nepoznanica, \mathbf{b} \in \mathbb{R}^{N} vektor desne strane čija je i-ta komponenta oblika \int_{0}^{1} f(x)x^{i}(1-x)dx te \mathbf{K} matrica sustava čiji (i,j) element ima oblik
\begin{align*} \int_{0}^{1} \left[\left(x^{j}(1-x)\right)'\left(x^{i}(1-x)\right)' + x^{j}(1-x)x^{i}(1-x)\right]dx &= \frac{ij}{i+j-1} - \frac{i(j+1) j(i+1)}{i+j} + \frac{(i+1)(j+1)}{i+j+1}\\\ &+ \frac{1}{i+j+1} - \frac{2}{i+j+2} + \frac{1}{i+j+3}. \end{align*}
Valja odmah primijetiti da je matrica koeficijenata \mathbf{K} popunjena elementima različitim od nule. No to nije sve, matrica je i loše uvjetovana što upućuje na ozbiljne opasnosti prilikom numeričkog rješavanja dobivenog sustava. U sljedećoj tablici je prikazana uvjetovanost matrice sustava za neke dimenzije potprostora.

\begin{array}{c|c}N& Cond(K)\\\ 3& 1.90 \cdot 10^{2} \\\ 5& 7.88 \cdot 10^{4}\\\ 7& 4.43 \cdot 10^{7}\\\ 10& 7.87 \cdot 10^{11}\\\ \end{array}

Kao posljednji argument koji ide u prilog numeričkoj nestabilnosti implementacije Galerkinove metode pogledajmo sljedeće slike.


 
 
Slika 1: Rezultati numeričke implementacije Galerkinove metode za f(x) = -x^{3} čije egzaktno rješenje je dano s u(x)=\frac{-7e^{x}}{e^{-1} - e}+\frac{7e^{-x}}{e^{-1} - e} - x^{3} - 6x.
a) Za n = 10, apsolutna pogreška je 7.87e\cdot 10^{11}, a uvjetovanost matice sustava je 4.02\cdot 10^{9}
b) Za n = 50, apsolutna pogreška je 1.11\cdot 10^{13} te uvjetovanost matice sustava je 5.28\cdot 10^{15}
c) Egzaktno rješenje
d) Popunjenost matrice krutosti za n = 25. Plava točkica označava ne-nul element u matrici na odgovarajućem mjestu.

 


Iako je s teorijske strane ovaj pristup jednostavan i elegantan, neke detalje ćemo ipak morati prilagoditi kako bismo osigurali spomenuta dobra svojstva i praktičnoj izvedbi.

3Metoda konačnih elemenata

U prethodnom primjeru prilikom numeričke izvedbe Galerkinove metode vidjeli smo da se rješavanje linearnog rubnog problema svodi na rješavanje sustava linearnih jednadžbi. No, matrica sustava je bila gusto popunjena i loše uvjetovana. Htjeli bismo prilagoditi Galerkinovu metodu tako da u konačnici dobijemo sustav koji je dobro uvjetovan i još po mogućnosti rijetko popunjen. Ipak, za računanje svakog elementa matrice potrebno je riješiti ponekad netrivijalan integral - što je vremenski skupo ukoliko je matrica gusto popunjena. Kada radimo sa sustavima čije su matrice rijetko popunjene otvara se još jedan bonus, mogućnost uporabe iterativnih metoda. Jedan način na koji možemo riješiti prethodna dva problema jest pažljivim odabirom baze konačnodimenzionalnog prostora V_{N}. Zathjevat ćemo da je broj baznih funkcija čiji se nosači sijeku s nosačem proizvoljne bazne funkcije što manji mogući. Odlučit ćemo se za prostor V_{N} koji se sastoji od funkcija koje su po dijelovima polinomi povezani s diskretizacijom domene \Omega. Jednostavno govoreći, Galerkinova metoda uz prethodne pretpostavke na bazne funkcije postaje metoda konačnih elemenata.
Konvergencija metode konačnih elemenata može se postići povećanjem stupnja polinoma u baznim funkcijama, profinjavanjem diskretizacije domene ili čineći obje istovremeno. Ilustrirajmo sada na primjeru navedene ideje.

Primjer 4. Promatramo
(6)
\begin{cases}-u'' + u = f \text{ na } \langle0,1\rangle\\ u(0) = 0, u'(1) = b \end{cases}
gdje je f\in L(\langle0,1\rangle) te b\in \mathbb{R}. Kao prije množenje s test funkcijom v \in V = \lbrace v \in H^{1}(\langle0,1\rangle) : v(0) = 0\rbrace, integriranjem po cijeloj domeni te primjenom parcijalne integracije daje slabu formulaciju
u\in V, \hspace{5mm} \int_{0}^{1}(u'v' + uv)dx = \int_{0}^{1}fvdx + bv(1), \hspace{5mm} \forall v\in V.
Opet, primjenom Lax-Milgramove leme slijedi da slaba formulacija ima jedinstveno rješenje.
Za N \in \mathbb{N} načinimo diskretizaciju domene I = [0, 1] tako da je 0 = x_{0} \lt x_{1} \lt ... \lt x_{N} = 1. Točke x_{i} nazivamo čvorovi, a podintervale I_{i} = [x_{i-1}, x_{i}] elementi. Označimo s h_{i} = x_{i} - x_{i-1} te s h = \max_{i = 1,...,N}h_{i}. Rješenje ćemo tražiti u prostoru
V_{h} = \lbrace v_{h} \in V : v_{h}|_{I_{i}}\in \mathcal{P}_{1}(I_{i}), i = 1,...,N\rbrace .
Zato što su po dijelovima glatke funkcije u H^{1}(I) ako i samo ako su u C(\overline{I}), ekvivaletno možemo definirati
V_{h} = \lbrace v_{h} \in C(\overline{I}) : v_{h}|_{I_{i}}\in \mathcal{P}_{1}(I_{i}), i = 1,...,N, v(0) = 0\rbrace .
Za bazne funkcije uvodimo tzv. hat funkcije
\begin{align*} \phi_{i}(x) = \begin{cases} (x - x_{i-1})/h_{i}, \hspace{5mm} &x_{i-1}\leq x \leq x_{i},\\ (x_{i+1}-x)/h_{i+1}, \hspace{5mm} &x_{i}\leq x \leq x_{i+1}, \\ 0, \hspace{5mm} &\text{inače}, \end{cases} \end{align*}
i posebno za desni rub kada je i = N,
\begin{align*} \phi_{i}(x) = \begin{cases} (x-x_{N-1})/h_{N}, \hspace{5mm} &x_{N-1}\leq x \leq x_{N}, \\ 0, \hspace{5mm} &\text{inače}, \end{cases} \end{align*}
Prethodne funkcije su neprekidne, po dijelovima glatke i očito linearno nezavisne. Lako se vidi da je prva slaba derivacija tako definiranih funkcija \phi_{i} po dijelovima konstantna funkcija.
 
Slika 2: Grafovi funkcija baznih funkcija


Neka je zadan prostor konačnih elemenata
V_{h} = span\lbrace \phi_{i} : i = 1,...,N\rbrace
uz problem
(7)
u_{h}\in V_{h} \hspace{5mm} \int_{0}^{1}(u_{h}'v_{h}' + u_{h}v_{h})dx = \int_{0}^{1}fv_{h}dx + bv_{h}(1) \hspace{5mm} \forall v\in V_{h},
koji prema Lax-Milgramovoj lemi ima jedinstveno rješenje. Kao i u prethodnom primjeru, ako zapišemo aproksimaciju rješenja kao linearnu kombinaciju elemenata baze,
u_{h} = \sum_{j = 1}^{N}u_{j}\phi_{j}
vidimo da se (7) može zapisati kao
\sum_{j=1}^{N}u_{j}\int_{0}^{1}(\phi_{i}'\phi_{j}' + \phi_{i}\phi_{j})dx = \int_{0}^{1}f\phi_{i}dx+b\phi_{i}(1), \hspace{5mm} i = 1,..., N.
Lako je provjeriti istinitost sljedećih formula, a njihovo uvažavanje će biti ključno za sastavljanje matrice sustava.
\begin{align*} \int_{0}^{1}\phi_{i}'\phi_{i-1}'dx &= -\frac{1}{h}, \hspace{5mm} i = 2,...,N,\\\ \int_{0}^{1}(\phi_{i}')^{2}dx &= \frac{2}{h}, \hspace{5mm} i =2,...,N-1,\\\ \int_{0}^{1}\phi_{i}\phi_{i-1}dx &= -\frac{h}{6}, \hspace{5mm} i = 2,...,N,\\\ \int_{0}^{1}(\phi_{i})^{2}dx &= \frac{2h}{3}, \hspace{5mm} i = 1,...,N-1,\\\ \int_{0}^{1}(\phi_{N}')^{2}dx &= -\frac{1}{h},\\\ \int_{0}^{1}(\phi_{N})^{2}dx &= -\frac{h}{3}. \end{align*}
Sada možemo problem (7) prikazati kao sustav linearnih jednadžbi \mathbf{Ku = b}.
\begin{align*} \mathbf{K} = \begin{bmatrix} \frac{2h}{3} + \frac{2}{h}& \frac{h}{6}- \frac{1}{h}&& \\\ \frac{h}{6}- \frac{1}{h}& \frac{2h}{3} + \frac{2}{h}& \frac{h}{6}- \frac{1}{h}& \\\ \ddots& \ddots& \ddots&& \\\& \frac{h}{6}- \frac{1}{h}& \frac{2h}{3} + \frac{2}{h}& \frac{h}{6}- \frac{1}{h} \\\&& \frac{h}{6}- \frac{1}{h}& \frac{h}{3} + \frac{}{h} \end{bmatrix} \end{align*}
je matrica sustava,
\mathbf{b} = \left(\int_{0}^{1}f\phi_{1}dx, \int_{0}^{1}f\phi_{2}dx,...,\int_{0}^{1}f\phi_{N-1}dx,\int_{0}^{1}f\phi_{N}dx + b\right)^{\tau}
vektor desne strane te
\mathbf{u} = (u_{1},...,u_{N})^{\tau}
vektor nepoznanica. Zahtijevali smo da bazne funkcije imaju što manji mogući nosač, posljedica je rijetko popunjena matrica sustava \mathbf{K} matrica. Dobra uvjetovanost matrice \mathbf{K} pokazana u sljedećoj tablici.



\begin{array}{ l | r } N& Cond(K)\\\ 3& 22.3962\\\ 5& 22.3962\\\ 7& 22.3962\\\ 10& 43.1613\\\ \end{array}

Na kraju primjera pogledajmo rezultate programske implementacije metode konačnih elemenata.


   
   
Slika 3: Rezultati numeričke implementacije metode konačnih elemenata za f(x) = -x^{3} uz rubne uvjete u(0) = u'(1) = 0 čije egzaktno rješenje je dano s u(x)=\frac{9e^{x}}{e^{-1} + e}-\frac{9e^{-x}}{e^{-1} + e} - x^{3} - 6x.
a) Egzaktno rješenje
b) Za n = 4, apsolutna pogreška je 0.0015, a uvjetovanost matice sustava je 32.64
c) Za n = 10, apsolutna pogreška je 9.62 \cdot 10^{-5} te uvjetovanost matice sustava je 150.25.
d) Popunjenost matrice \mathbf{K} za n = 25. Plava točkica označava ne-nul element na odgovarajućem mjestu u matrici.

4Još jedan primjer

U nastavku promotrimo jednadžbu koja opisuje otklon žice klavira pod utjecajem sile f. Žica klavira ima dva dijela; žičanu jezgru te žicu namotanu oko nje. Uzevši to u obzir, jasno je da se takav objekt ne ponaša isključivo kao tanka žica. Žica klavira vibrira dijelom zbog napetosti, a dijelom zbog same čvrstoće žice. Iz tog razloga modelu dodajemo član koji opisuje ponašanje tankog štapa.

(8)
\begin{cases}u^{(4)} - u^{''} = f \text{ na } \langle0,1\rangle\\ u(0) = u'(0) = u(1) = u'(1) = 0. \end{cases}

Slaba formulacija problema je

u\in V, \hspace{5mm} \int_{0}^{1}(u''v'' + u'v')dx = \int_{0}^{1}fvdx \hspace{5mm} \forall v\in V,

gdje je V = H_{0}^{2}(\langle0,1\rangle). Ako želimo odabrati V_{h}, takav da je V_{h} \subset V, onda svaka funkcija u V_{h} mora bit barem klase \mathcal{C}^{1}. Pretpostavimo da se V_{h} sastoji od po dijelovima polinoma stupnja manjeg ili jednakog p. Uvjet da su funkcije u V_{h} klase \mathcal{C}^{1} nameće 2(N-1) uvjeta, a rubni uvjeti još 4 uvjeta. Sve skupa imamo

dim(V_{h}) = (p+1)N - 2(N-1) - 4 = (p-1)N - 2

stupnja slobode da odaberemo bazne funkcije sa što manjim nosačem. Jasno je da je najmanji dozvoljeni p = 2. Međutim odlučit ćemo se za p = 3 da bismo povećali broj stupnjeva slobode za odabir baznih funkcija. Za p =3 prostor konačnih elemenata je

V_{h} = \lbrace v_{h}\in \mathcal{C}^{1}(\overline{I}) : v_{I_{i}}\in \mathcal{P}_{3}(I_{i}), I =1,...,N, v_{h}(x) = v_{h}'(x) = 0 \text{ za } x = 0,1\rbrace .

Možemo konstruirati bazne funkcije s malim nosačima koristeći uvjete interpolacije

\begin{align*} \phi_{i}(x_{j}) = \delta_{ij}, \hspace{4mm} \phi_{i}'(x_{j}) = 0,\\ \psi_{i}(x_{j}) = 0, \hspace{4mm} \psi_{i}'(x_{j}) = \delta_{ij}. \end{align*}

Kako bismo izbjegli tehničke komplikacije, poželjno je raditi na referentnom elementu. Ideja je sav posao obaviti na jednom elementu, tzv. referentnom elementu, a onda postignuto preslikati na ostale elemente. Neka je I_{0} = \langle 0, 1\rangle, preslikavanje

F_{i} : I_{0} \to I_{i}, \hspace{5mm} F_{i}(\xi) = x_{i-1}+h_{i}\xi

je bijekcija. Na referentnom elementu I_{0} konstruiramo polinome trećeg stupnja koji zadovoljavaju interpolacijske uvjete

\begin{align*} \Phi_{0}(0) = 1, \hspace{4mm} \Phi_{0}(1) = 0,\hspace{4mm} \Phi_{0}'(0) = 0,\hspace{4mm} \Phi_{0}'(1) = 0,\\\ \Phi_{1}(0) = 0, \hspace{4mm} \Phi_{1}(1) = 1,\hspace{4mm} \Phi_{1}'(0) = 0,\hspace{4mm} \Phi_{1}'(1) = 0,\\\ \Psi_{0}(0) = 0, \hspace{4mm} \Psi_{0}(1) = 0,\hspace{4mm} \Psi_{0}'(0) = 1,\hspace{4mm} \Psi_{0}'(1) = 0,\\\ \Psi_{1}(0) = 0, \hspace{4mm} \Psi_{1}(1) = 0,\hspace{4mm} \Psi_{1}'(0) = 0,\hspace{4mm} \Psi_{1}'(1) = 1. \end{align*}

Jednostavnim računom određujemo da su tražene funkcije

\begin{align*} &\Phi_{0} (\xi) = (1+2\xi)(1-\xi)^{2},\\ &\Phi_{1} (\xi) = (3-2\xi)\xi^{2}, \\ &\Psi_{0}(\xi) = \xi(1-\xi)^{2},\\ &\Psi_{1}(\xi) = -(1-\xi)\xi^{2}. \end{align*}

Sad možemo odrediti bazne funkcije

\begin{align*} \phi_{i} = \begin{cases} \Phi_{1}(F_{i}^{-1}(x)), \hspace{5mm} &x\in I_{i},\\ \Phi_{0}(F_{i+1}^{-1}(x)), &x\in I_{i+1},\\ 0, &\text{inače}, \end{cases} \end{align*}

i

\begin{align*} \psi_{i} = \begin{cases} h_{i}\Psi_{1}(F_{i}^{-1}(x)), \hspace{5mm} &x\in I_{i},\\ h_{i+1}\Psi_{0}(F_{i+1}^{-1}(x)), &x\in I_{i+1},\\ 0, &\text{inače}, \end{cases} \end{align*}

Kao i u prethodnim primjerima, želimo formulirati sustav linearnih jednadžbi. Koeficijente matrice sustava računat ćemo preko referentnog elementa, npr.

\begin{align*} k_{i-1,i} = \int_{0}^{1}(\phi_{i-1})''(\phi_{i})''+ (\phi_{i-1})'(\phi_{i})'dx = \int_{I_{i}}(\phi_{i-1})''(\phi_{i})'' + (\phi_{i-1})'(\phi_{i})'dx \end{align*}

koristeći definiciju baznih funkcija i funkcije F_{i} možemo izračunati

\begin{align*} k_{i-1,i} &= \frac{1}{h_{i}^{3}}\int_{I_{0}}6(2\xi-1)6(1-2\xi)d\xi + \frac{1}{h_{i}}\int_{I_{0}}6\xi(\xi-1)6\xi(1-\xi)d\xi\\ &= -\frac{12}{h_{i}^{3}} - \frac{6}{5h_{i}}. \end{align*}

Ostale elemente matrice možemo izračunati analogno, a tako dobivena matrica će biti petdijagonalna. Komentirajmo za kraj da je za višedimenzionalne probleme tehnika referentnog elementa neophodna, kako za teorijsku analizu greške tako i za praktičnu implementaciju. Valja primijetiti da smo u prethodnom jednodimenzionalnom primjeru zahtijevali da funkcije iz V_{h} budu polinomi trećeg stupnja kako bismo postigli da funkcije iz V_{h} \subset V budu klase \mathcal{C}^{1}. U višim dimenzijama je potreban još veći stupanj polinoma kako bismo to postigli. Jedan način kako je moguće riješiti taj problem jest korištenjem nekonformnih elemenata ([5]).
Za sam kraj, ostavljamo čitatelju jedan programerski podhvat.


Zadatak. Napravite programsku izvedbu metode konačnih elemenata za sljedeći problem

(9)
\begin{cases}u^{(4)} - u^{''} = -x^{3} \text{ na } \langle0,1\rangle\\ u(0) = u'(0) = u(1) = u'(1) = 0. \end{cases}

Egzaktno rješenje ovog problema je

\begin{align*}u(x) =& \frac{1}{20(e-3)(e-1)}\cdot e^{-x}(e^{x+2}(x^{5}+20x^{3}-44x+44) - 2e^{x+1}(2x^{5}+40x^{3}-65x+44)\\\ &+e^{x}(3x^{5}+60x^{3}-86x-86)+86e^{2x}-21e^{2x+1}-44e^{2}+109e). \end{align*}

Zahvaljujem se prof.dr.sc. Josipu Tambači na detaljnom čitanju ovog članka i korisnim savjetima.

5Literatura

[1] Atkinson, Han. Theoretical Numerical Analysis: a functional framework. Vol 19. Springer, 2009.
[2] Dhatt, Gouri i Touzot, Gilbert. Finite Element Method. John Willey and Sons, 2012.
[3] Brenner, Susanne i Carstensen, Carsten. Encyclopedia of Computational Mechanics. Chapter 4, Finite Element Methods. John Wiley and Sons, 2004.
[4] Jurak, Mladen. Praktikum primijenjene matematike II. Skripta, PMF Matematički odsjek, 2008.

[5] Ciarlet, Philippe G. The finite element method for elliptic problems. Elsevier, 1978.

 

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.)
 

 

Stabilnost rotirajućeg štapa

 
Josip Tambača
PMF-Matematički odjel
Sveučilište u Zagrebu
tambaca@math.hr
Josip Tokmačić
tokma80@net.hr


Autori zahvaljuju prof. dr. sc. Branimiru Matijaševiću i prof. dr. sc. Zvonimiru Guzoviću na organizaciji i izvođenju eksperimenta čiji su rezultati navedeni u posljednjem odjeljku rada.


1Uvod

Tijela koja rotiraju nalazimo u mnogim mehaničkim uređajima, pa i onima kojima se i sami služimo. Ovdje ćemo se posebno baviti izduženim tijelima koja se rotiraju oko svoje „izdužene” osi. Primjer su svrdla raznih bušilica (klasične ručne ili one za bušenja u geologiji, vađenju nafte), razne osovine (vratila) za prijenos snage, kao kod automobila, brodova ili u primjeru rotora turbina elektrana.


Slika 1: Osovina brodskog vijka


Slika 2: Razne bušilice, ne samo ručne


Slika 3: Osovine koje spajaju generatore i lopatice turbina


Slika 4: Rotor parne turbine


Kod svih ovih primjera mehaničkih sustava tijelo rotira oko svoje osi. Potencijalni problemi dolaze zbog utjecaja centrifugalne sile koja teži izbaciti sustav iz pravilne vrtnje.


Slika 5: Deformirani položaj štapa

Tipično, to se događa za dovoljno velike brzine vrtnje. Posljedice su velike vibracije koje mogu uništiti čitav sustav. Jednom kad izgradimo konkretan sustav, relativno ga je lako eksperimentalno testirati i provjeriti radi li on korektno u radnim uvjetima. No, ako smo ga izgradili „na slijepo”, velike su šanse da nije racionalno izgrađen (ili je predimenzioniran ili iskazuje određene nedostatke, primjerice nestabilnost). Stoga je od interesa modelirati konkretni mehanički sustav da bismo ga lakše projektirali. U našem konkretnom slučaju poželjno je unaprijed znati koje su to brzine koje će producirati nestabilnu vrtnju.

U ovom nam je radu u cilju izložiti čitav proces modeliranja rotirajućih štapova. Počinjemo s opisom teorije i smještanjem problema u okvire te teorije. Zatim tu teoriju primjenjujemo na konkretni problem rotirajućeg štapa. Dobiveni model smještamo u matematički kontekst, te matematički analiziramo problem. Nakon toga određujemo analitička rješenja problema za neke specijalne slučajeve (i jedine kada je to i moguće u našem primjeru), te donosimo numeričke aproksimacije rješenja. Na kraju uspoređujemo dobivene rezultate (analitičke i numeričke) s rezultatima eksperimenta, te time konačno potvrđujemo model.

2Model elastičnog štapa

2.1Geometrija štapa

„Izduženo” tijelo koje se rotira u nastavku nazivamo tanki štap. Tanki štap \Omega \subset \mathbb{R}^{3} opisujemo kao trodimenzionalno tijelo kojem su dvije dimenzije („debljine”) malene u odnosu na treću („duljinu”, L\gt 0):


Slika 6: Tanki ravni štap

u inženjerskoj literaturi tipično se pretpostavlja da je debljina / duljina \lt 1/ 30. Malo preciznije, neka je dana familija \lbrace D(x_{3})\subset \mathbb{R}^{2}_{x_{1} x_{2}}: x_{3} \in (0, L)\rbrace. Štap definiramo s
(1)
\Omega = \left\lbrace x=(x_{1}, x_{2}, x_{3}) \in \mathbb{R}^{3} : (x_{1}, x_{2}) \in D(x_{3}), x_{3} \in (0, L)\right\rbrace .
Skupove D(x_{3}) zovemo poprečni presjeci štapa u točki x_{3}. Dopuštamo ovisnost o točki poprečnog presjeka jer su nam baš takvi i realni primjeri. Sada tankoću štapa iskazujemo zahtjevom: \mathop{\rm diam}\nolimits{D(x_{3})}\ll L. Specijalno, od interesa će nam biti samo ravni štap. To je onaj kojemu je centralna linija (linija koja prolazi kroz težišta poprečnih presjeka) ravna. Stoga pretpostavljamo da vrijedi
\mathop{\int\!\!\!\int}\limits_{D(x_{3})} x_{\alpha} \, dx_{1} dx_{2} =0, \quad \alpha=1, 2, \quad x_{3} \in (0, L).
Radi jednostavnosti, ovdje dodatno pretpostavljamo da poprečni presjeci zadovoljavaju dodatni uvjet simetrije
\mathop{\int\!\!\!\int}\limits_{D(x_{3})} x_{1} x_{2} \, dx_{1} dx_{2} = 0, \qquad x_{3} \in (0, L)
(zadovoljen je, naprimjer, za kružni i pravokutni poprečni presjek).

2.2Elastičnost

Djelujući silom naprimjer na tanki čelični štap, on se deformira. Uklonimo li silu kojom smo djelovali, štap se vraća u nedeformirani (prvobitni) položaj. Tipičan primjer takvog ponašanja je opruga. Tijela koja se pri djelovanju sile deformiraju, a nakon uklanjanja sile vraćaju u prvobitni položaj, nazivamo elastičnim tijelima, a teoriju koja se njima bavi teorijom elastičnosti.

Tanki elastični štap je prema prethodnom opisu trodimenzionalno elastično tijelo. Jednadžbe trodimenzionalne teorije elastičnosti, čak i one linearizirane (pojednostavljene uz pretpostavku neznatne deformacije), komplicirane su. No, u slučaju posebnih geometrija tijela, kao u slučaju štapa, moguće je izvesti jednostavnije modele ponašanja (slično je u primjeru opruge, koja je trodimenzionalno elastično tijelo; njeno je ponašanje vrlo kompleksno za precizan opis, ali znamo da je „osnovno” ponašanje modelirano jednadžbom F=k x; F je sila koju primjenimo na oprugu, x pripadno izduženje (skraćenje) opruge, a k krutost opruge). Geometrijski tanki štap možemo „dobro” opisati njegovom centralnom krivuljom, parametriziranom s (0,0,x_{3}), x_{3}\in (0,L), a „dobro” ovdje znači do na najveći dijametar poprečnih presjeka. Pokazuje se da se i deformacija tankog elastičnog štapa može opisati deformacijom centralne krivulje štapa, parametriziranom s \mathbf{\Psi}(x_{3}) \in \mathbb{R}^{3}, x_{3}\in (0,L).

I opet je pogreška takvog modela, među ostalim, vezana za najveći dijametar poprečnih presjeka. Izvod takvog modela daje i malo više. Naime, pokazuje se da se deformacija tankog elastičnog štapa može opisati kroz deformaciju centralne linije štapa, te kroz infinitezimalnu krutu deformaciju poprečnih presjeka (infinitezimalna kruta deformacija je translacija i rotacija, s tim da je rotacija aproksimirana sumom matrica I+A, gdje je I identiteta, a A antisimetrična matrica. Naime, može se pokazati da je svaka rotacija R=e^{A}, za neku antisimetričnu matricu, pa je I+A zapravo Taylorov polinom 1. stupnja za R).


Slika 7: Pomak tankog štapa

Takva deformacija \mathbf{\Psi}(x), x \in \Omega s pomoću pomaka \mathbf{U}(x) = \mathbf{\Psi}(x) - x može se zapisati formulama
(2)
\begin{eqnarray} U_{1} (x) &=& u_{1} (x_{3}) - x_{2} \omega (x_{3}),\\ U_{2} (x) &=& u_{2} (x_{3}) + x_{1} \omega (x_{3}),\\ U_{3} (x) &=& u_{3} (x_{3}) - x_{1} u'_{1} (x_{3}) - x_{2} u'_{2} (x_{3}). \end{eqnarray}
Ovdje je \mathbf{u}(x_{3})=(u_{1}(x_{3}),u_{2}(x_{3}),u_{3}(x_{3})) pomak točke (0,0,x_{3}) centralne linije štapa, u'_{1}(x_{3}) opisuje kut zakreta poprečnog presjeka na mjestu x_{3} oko osi \mathbf{e}_{2}, u'_{2}(x_{3}) opisuje kut zakreta poprečnog presjeka na mjestu x_{3} oko osi \mathbf{e}_{1}, a \omega(x_{3}) opisuje kut zakreta poprečnog presjeka na mjestu x_{3} oko osi \mathbf{e}_{3} (torziju). Ovakav oblik pomaka naziva se Bernoulli–Navierov (ili Kirchhoff–Love) pomak.

Budući da ćemo promatrati stabilnost rotirajućeg štapa u odnosu na transverzalni pomak (okomit na centralnu liniju), u nastavku pišemo samo jednadžbe koje zadovoljavaju transverzalni pomaci (u_{1},u_{2}). To su obične diferencijalne jednadžbe četvrtog reda:
(3)
\begin{eqnarray} &&E \left(I_{1} u''_{1}\right)'' - \left(\sigma u_{1}'\right)' - f_{1} + m_{2}' =0,\\ &&E \left(I_{2} u''_{2}\right)'' - \left(\sigma u_{2}'\right)' - f_{2} - m_{1}' =0. \end{eqnarray}
(vidi npr. [4], [6]). Ovdje su (f_{1}, f_{2}) i (m_{1}, m_{2}), redom, linijska gustoća transverzalne sile i momenta (to znači da je \int_{a}^{b} f_{1} dx_{3} komponenta ukupne sile u smjeru \mathbf{e}_{1}=(1,0,0) na komad štapa \lbrace x\in \mathbb{R}^{3}: (x_{1},x_{2} \in D(x_{3})), x_{3} \in [a,b]\rbrace). E\gt 0 je Youngov modul elastičnosti materijala (opisuje svojstvo materijala od kojeg je štap napravljen), a I_{\alpha} (x_{3}), \alpha=1, 2, su momenti inercije poprečnog presjeka D(x_{3}) – opisuju geometriju poprečnog presjeka štapa:
I_{\alpha} (x_{3}) = \mathop{\int\!\!\!\int}\limits_{D(x_{3})} x_{\alpha}^{2} \, dx_{1} dx_{2}, \quad \alpha = 1, 2.
Funkcija \sigma(x_{3}) je longitudinalna napetost štapa (npr. u slučaju bušilice sila kojom je pritišćemo primjerice na zid). Strogi izvod ovog modela može se naći u [6] i [7].

3Rotirajući elastični štap

U nastavku ćemo iskoristiti model (3) za štap koji se rotira konstantnom kutnom brzinom \gamma, te na njega ne djeluju druge sile osim centrifugalne. Očito je riječ o dinamičkom problemu, što znači da bismo se morali koristiti jednadžbama gibanja, a ne ravnoteže. No, ako odaberemo koordinatni sustav vezan za štap koji se rotira, u tom koordinatnom sustavu štap miruje. U modelu (3) transverzalni pomak uzrokovan je centrifugalnom silom i njome uzrokovanim momentom. Stoga ih, u koordinatnom sustavu u kojem štap miruje, računamo po formulama
(4a)
\mathbf{f}(x_{3}) = - \varrho \mathop{\int\!\!\!\int}\limits_{D(x_{3})} \gamma \mathbf{e}_{3} \times \left(\gamma \mathbf{e}_{3} \times \left( x + \mathbf{U}(x)\right)\right) \, dx_{1} dx_{2},
(4b)
\mathbf{m}(x_{3}) = - \varrho \mathop{\int\!\!\!\int}\limits_{D(x_{3})} \left(x + \mathbf{U}(x) - (x_{3} \mathbf{e}_{3} + \mathbf{u}(x_{3}))\right) \times\left(\gamma \mathbf{e}_{3} \times \left(\gamma \mathbf{e}_{3} \times \left(x + \mathbf{U}(x)\right)\right)\right) \, dx_{1} dx_{2},
gdje je \mathbf{e}_{3} =(0, 0, 1); \varrho \in \mathbb{R}^{+} je linijska gustoća mase (\rho d je masa komada štapa duljine d), a \times označava vektorski produkt (ovo su uobičajene formule za centrifugalnu silu \mathbf{\gamma} \times (\mathbf{\gamma} \times \mathbf{r}), gdje je \mathbf{\gamma} vektor kutne brzine, a \mathbf{r} radij vektor mase; jedini je dodatak integral po svim točkama tijela jer nije riječ o masi u točki).

Koristeći Bernoulli-Navierov oblik pomaka (2), u stanju smo prikazati silu i moment s pomoću transverzalnog pomaka centralne linije:
(5)
f_{1} = \varrho S \gamma^{2} u_{1}, \quad f_{2} = \varrho S \gamma^{2} u_{2},\quad m_{1} = \varrho I_{2} \gamma^{2} u_{2}', \quad m_{2} = - \varrho I_{1} \gamma^{2} u_{1}';
ovdje je S(x_{3}) = | D(x_{3})| površina poprečnog presjeka na mjestu x_{3}.

Uvrštavanjem ovih sila i momenata u jednadžbe modela (3) dobivamo dvije jednadžbe
\left(I_{\alpha} u''_{\alpha}\right)'' - \frac{1}{E} \left(\sigma u'_{\alpha}\right)' - \lambda \left(\left( I_{\alpha} u'_{\alpha}\right)'+S u_{\alpha} \right) = 0,\qquad \alpha=1,2,
gdje smo uveli konstantu
(6)
\lambda = \frac{\varrho \gamma^{2}}{E}.
Ovo su jednadžbe za transverzalni pomak u oba transverzalna smjera. S obzirom na to da su jednadžbe iste, u nastavku pišemo samo jednu
(7)
\left(I_{\alpha} u''\right)'' - \frac{1}{E} \left(\sigma u'\right)' - \lambda \left(\left( I u'\right)'+S u\right) = 0.
U slučaju kada su oba kraja štapa ukliještena (nije dopušten progib ni rotacija poprečnog presjeka na rubu), uz ovu jednadžbu pridružujemo i rubne uvjete
(8)
u = u' =0 \text{ za } x=0, L.


Slika 8: Oblik deformacije štapa za ukliješteni štap


U slučaju kad su oba kraja učvršćena šarnirom (nije dopušten progib, dok je rotacija slobodna na krajevima), rubni uvjeti su
(9)
u = u'' =0 \text{ za } x=0, L.


Slika 9: Oblik deformacije štapa za štap učvršćen šarnirom

4Kritična brzina štapa

Zadaća za pomak štapa pri rotaciji sada glasi: naći funkciju u koja zadovoljava jednadžbu (7) te jedan od rubnih uvjeta (8) ili (9), ovisno o problemu koji promatramo (čak možemo i kombinirati različite rubne uvjete na oba kraja). Jasno je da za svaki \lambda \in \mathbb{R} jedno rješenje ovih zadaća funkcija u=0. Pitanje je postoje li možda i neka druga rješenja za određene \lambda (možda ne za sve). Stoga je i \lambda nepoznanica u problemu i on sada glasi: naći sve \lambda \in \mathbb{R} i funkcije u koje nisu nul-funkcije za koje vrijedi
(10)
\begin{split} &\left(I u''\right)'' - \frac{1}{E} \left(\sigma u'\right)' = \lambda \left(\left( I u'\right)'+S u\right),\\ & u(0) = u'(0)= u(L) = u'(L) =0, \end{split}
ili u slučaju drugih rubnih uvjeta
(11)
\begin{split} &\left(I u''\right)'' - \frac{1}{E} \left(\sigma u'\right)' = \lambda \left(\left( I u'\right)'+S u\right),\\ & u(0) = u''(0)= u(L) = u''(L) =0. \end{split}
Zadaće ovog tipa nazivaju se svojstvene zadaće (ili spektralne zadaće), a ovo su primjeri svojstvenih zadaća za linearni diferencijalni operator (4. reda). Rješenja su zapravo uređeni parovi (\lambda, u_{\lambda}) realnog broja \lambda i nenul-funkcije u_{\lambda} koja rješava rubnu zadaću za diferencijalnu jednadžbu za dani \lambda. \lambda se naziva svojstvena vrijednost, a pripadni u_{\lambda} svojstveni vektor ili svojstvena funkcija. Za primijetiti je da takav u_{\lambda} nije jedinstven jer je i svaka funkcija oblika c u_{\lambda}, za c\in \mathbb{R} također svojstvena funkcija.

Najjednostavniji primjer svojstvenih zadaća su svojstvene zadaće za linearne operatore na konačnodimenzionalnim vektorskim prostorima koje se reduciraju na matrične svojstvene zadaće (vidi Linearnu algebru).

Matematička teorija ovakvih svojstvenih zadaća bazirana je na spektralnoj teoriji kompaktnih linearnih operatora [3]. Osnovni je rezultat za našu svojstvenu zadaću da je spektar (skup svih svojstvenih vrijednosti) prebrojiv (ima članova koliko ima i prirodnih brojeva), te da nema konačnog gomilišta. Ovo zajedno daje da je spektar diskretan. Spektar u nastavku označavamo s \Lambda.

U ovom problemu nenul–rješenja dolaze iz dvaju izvora: rotacije štapa te longitudinalne sile koja djeluje na štap. Naprimjer, pretpostavimo da se štap ne rotira; uzmemo li dovoljno tanki čelični štap i pritisnemo li ga duž osi štapa, doći će do transverzalnog progiba ako smo primijenili dovoljno veliku silu.

Taj slučaj nestabilnosti bez vrtnje (Eulerova nestabilnost, vidi [4]) ovdje isključujemo zahtjevom (vidi [1])
\sigma_{\min{}} + E I_{\min{}} \left( \frac{\pi}{L}\right)^{2} \gt 0,
pri čemu je
\sigma_{\min{}} = \min_{x_{3} \in [0, L]}{\sigma(x_{3})}, \quad I_{\min{}} ={\min_{x_{3} \in [0, L]}{I (x_{3})}}.

Zbog definicije broja \lambda u (6) i našeg problema jasno je da su nam od interesa samo \lambda \gt 0. Stoga se ograničavamo na promatranje skupa \Lambda \cap \mathbb{R}^{+}. Označimo
\lambda_{0} = \min \Lambda \cap \mathbb{R}^{+},
s tim da definiramo \lambda_{0} = +\infty ako je skup \Lambda \cap \mathbb{R}^{+} prazan. Sada je štap stabilan za kutne brzine rotacije \gamma koje zadovoljavaju
0 \lt \gamma \lt \gamma_{\scriptsize{\mbox{kritična}}} = \left( \frac{E\lambda_{0}}{\varrho}\right)^{1/2}.
Kutnu brzinu vrtnje \gamma_{\scriptsize{\mbox{kritična}}} zovemo kritična brzina. Osim u jednom, dva slučaja kritičnu brzinu nije moguće točno izračunati. Ostaju nam dvije mogućnosti aproksimacije kritične brzine. Jedna je bazirana na analitičkoj ocjeni kritične brzine odozdo i izvedena je u [1]. Druga je numerička aproksimacija bazirana na metodi konačnih elemenata (Finite Element Method, FEM) i prikazana je u [5]. Analitičku ocjenu odozdo donosimo sljedećem teoremu

Teorem 1.([1]) Neka je S_{\max} = \max_{x_{3} \in [0, L]}{S(x_{3})}. Tada vrijedi
\begin{eqnarray*} &&S_{\max} \leq \left( \frac{\pi}{L}\right)^{2} I_{\min},\qquad \Rightarrow \qquad \gamma_{\scriptsize{\mbox{crit}}} = \infty,\\ &&S_{\max} \gt \left( \frac{\pi}{L}\right)^{2} I_{\min},\qquad \Rightarrow \qquad \gamma_{\scriptsize{\mbox{crit}}} \geq \overline{\gamma}= \frac{\pi}{L} \varrho^{-1/2} \left( \frac{\sigma_{\min} + E I_{\min} \left( \frac{\pi}{L}\right)^{2}}{S_{\max} - I_{\min} \left( \frac{\pi}{L}\right)^{2}}\right)^{1/2}. \end{eqnarray*}

Ovaj rezultat zapravo kaže: ako je štap dovoljno „debeo” nema nestabilnosti; inače, nestabilnost je moguća samo za brzine rotacije veće od \overline{\gamma}.

5Primjeri

Ovdje donosimo dva primjera. U prvom primjeru (vjerojatno je i jedini takav) znamo izračunati analitička rješenja našeg problema. U drugom primjeru svedemo određivanje svojstvenih vrijednosti na problem nultočaka nelinearne funkcije. U oba slučaja riječ je o štapu konstantnog poprečnog presjeka te konstantne uzdužne sile \sigma. Za sve druge slučajeve zadaću moramo numerički rješavati.

Primjer 2.([2), str. 116.] Neka je poprečni presjek konstantan. Tada su I, S konstante. Pretpostavimo nadalje da je \sigma \geq 0 konstanta. Označimo \theta = \frac{\sigma}{E I}. Tada zadaća (11) poprima oblik
(12)
\begin{split} &u'''' - (\theta + \lambda ) u'' - \lambda \frac{S}{I} u = 0,\\ &u(0)=u''(0)=u(L)=u''(L)=0. \end{split}
Ovu diferencijalnu jednadžbu možemo riješiti eksplicitno (naći opće rješenje, tj. odrediti skup svih rješenja), te iskoristiti rubne uvjete da bismo riješili i svojstvenu zadaću (tj. skup nenul–rješenja zadaće). Više o diferencijalnim jednadžbama te rješavanju običnih diferencijalnih jednadžbi možete naći u [8]. Karakteristična jednadžba pripadna jednadžbi (12) glasi
k^{4} - (\theta + \lambda) k^{2} - \lambda \frac{S}{I} =0.
Za nultočke ove jednadžbe vrijedi
(k^{2})_{1/2} = \frac{1}{2} (\theta + \lambda) \pm \sqrt{\frac{1}{4} (\theta + \lambda)^{2} + \lambda \frac{S}{I}}.
Budući da tražimo nenegativne svojstvene vrijednosti, označimo li
(13)
\alpha = \left( \frac{1}{2} (\theta + \lambda) + \sqrt{\frac{1}{4} (\theta + \lambda)^{2} + \lambda \frac{S}{I}}\right)^{1/2},\quad \beta = \left( -\frac{1}{2} (\theta + \lambda) + \sqrt{\frac{1}{4} (\theta + \lambda)^{2} + \lambda \frac{S}{I}}\right)^{1/2},
rješenja karakteristične jednadžbe su \alpha, - \alpha, i \beta, - i \beta, pa je opće rješenje jednadžbe (12) oblika
u(z) = A \mathop{\rm sh}\nolimits{\alpha z} + B \mathop{\rm ch}\nolimits{\alpha z} + C \sin{\beta z} + D \cos{\beta z},
jer su \alpha \ne 0 i \beta \ne 0. Iz rubnih uvjeta (12) na lijevoj bazi štapa slijedi B=D=0, dok iz rubnih uvjeta na desnoj bazi slijedi
\begin{eqnarray*} 0 &=& A \mathop{\rm sh}\nolimits{\alpha L} + C \sin{\beta L},\\ 0 &=& A \alpha^{2} \mathop{\rm sh}\nolimits{\alpha L} - C \beta^{2} \sin{\beta L}. \end{eqnarray*}
Ovaj sustav ima netrivijalna rješenja ako je determinanta sustava jednaka nuli:
-(\alpha^{2} + \beta^{2}) \mathop{\rm sh}\nolimits{\alpha L} \sin{\beta L} =0.
Slijedi da svojstvenih vrijednosti \lambda ima prebrojivo (kao što smo i znali) te da su one određene jednadžbom
\beta_{n} L = n \pi, \quad n \in \mathbb{N}.
Stoga su svojstvene vrijednosti dane s
\lambda_{n} = \left( \frac{n \pi}{L}\right)^{2} \frac{\frac{\sigma}{E I}\left( \frac{n\pi}{L}\right)^{2}}{\frac{S}{I} - \left( \frac{n \pi}{L}\right)^{2}}, \quad n \in \mathbb{N}.
Najmanja svojstvena vrijednost dobiva se za n=1 te je kritična brzina dana izrazom
\gamma_{\scriptsize{\mbox{kritična}}} = \gamma_{1} = \sqrt{\frac{E \lambda_{1}}{\varrho}} = \frac{\pi}{L} \sqrt{\frac{\frac{\sigma}{I} + E \left( \frac{\pi}{L}\right)^{2}}{\varrho \left( \frac{S}{I} - \left( \frac{\pi}{L}\right)^{2}\right)}}.

Označimo
(14)
c = \frac{\sigma L^{2}}{E I}, \quad d = \frac{I}{S L^{2}}\lt \frac{1}{\pi^{2}}.
Sada kritičnu brzinu možemo zapisati s
(15)
\gamma_{\scriptsize{\mbox{kritična}}} = \frac{1}{L} \left( \frac{E}{\varrho}\right)^{1/2} \gamma_{0}, \qquad \gamma_{0} = \pi \left(d \frac{c + \pi^{2}}{1 - d \pi^{2}}\right)^{1/2}.
Graf od \gamma_{0} u ovisnosti o parametrima c i d dan je na Slici 10.


Slika 10: Slika grafa \gamma_{0} u ovisnosti o c i d

Parametar c možemo interpretirati kao uzdužnu silu i nalazi se na osi x, dok parametar d možemo interpretirati kao „debljinu” štapa i nalazi se na osi y (argument ide do 1/\pi^{2}). Zaključujemo da što je štap deblji, kritična brzina je veća, pa je štap stabilniji. Također, što je uzdužna sila veća, dobivamo stabilniji štap.

U Tablici 1 vrijednosti \gamma_{\scriptsize{\mbox{kritična}}} za čelični štap duljine L=1\,{\text{m}} kružnog poprečnog presjeka dane su u ovisnosti o radijusu štapa R. I ovdje vidimo da „deblji” štap postaje stabilniji. Brzinu \gamma koja odgovara 2. svojstvenoj vrijednosti nazivamo 2. kritična brzina, te tako redom dalje. I ove su brzine od interesa jer se ponekad radna brzina vrtnje sustava nalazi između 1. i 2. kritične brzine.


R[m] 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
\gamma_{\scriptsize{\mbox{kritična}}}[rpm] 2440 4883 7329 9781 12240 14708 17187 19679 22186 24711
Tablica 1: Kritične brzine u okretajima u minuti za različite radijuse štapa


Slika 11: Položaj pri 1. kritičnoj brzini


Slika 12: Položaj pri 2. kritičnoj brzini


Slike položaja pri prve dvije kritične brzine su na Slikama 11, 12.

Primjer 3.[[2], str. 116.] Neka je poprečni presjek konstantan. Tada su I i S konstante. Uz oznaku \theta = \frac{\sigma}{E I} i pretpostavku \sigma\geq 0 jednadžba (10) poprima oblik
(16)
\begin{split} &u'''' - (\theta + \lambda ) u'' - \lambda \frac{S}{I} u = 0,\\ &u(0)=u'(0)=u(L)=u'(L)=0. \end{split}
U ovom slučju ne možemo izračunati kritičnu brzinu analitički, već je određena najmanjom pozitivnom nultočkom eksplicitno dobivene nelinearne funkcije. Potpuni račun možete naći u [2].
Slike položaja pri prve dvije kritične brzine su na Slikama 13, 14. Obratite pozornost na ponašanje štapa na krajevima.


Slika 13: Položaj pri 1. kritičnoj brzini

U Tablici 2 vrijednosti \gamma_{\scriptsize{\mbox{kritična}}} za ravni čelični štap duljine L=1m s kružnim poprečnim presjekom dane su u ovisnosti o radijusu R. Uspoređujući ove vrijednosti s Tablicom 1, vidimo da je u ovom primjeru štap stabilniji, što je i očekivano.


R[m] 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
\gamma_{\scriptsize{\mbox{crit}}}[rpm] 5533 11071 16619 22183 27768 33378 39020 44699 50421 56192
Tablica 2: Kritične brzine u okretajima u minuti u odnosu na radijus R štapa


Slika 14: Položaj pri 2. kritičnoj brzini


6Numerički rezultati

Model rotirajućeg štapa nije moguće analitički riješiti čim neki od koeficijenata nije konstantan (npr. poprečni presjek nije konstantnog radijusa duž štapa). U svim tim slučajevima problem rješavamo numerički, metodom konačnih elemenata, vidi [5].

Primjer 4.[Varijabilni poprečni presjek] U ovom primjeru promatramo štap duljine L=1m, radijusa poprečnog presjeka
R(x) = \left\lbrace \begin{array}{ll} 0.01\text{m} & x\in [0, 0.4\rangle \cup \langle 0.6, 1],\\ \tilde{r} & x \in [0.4,0.6], \end{array} \right.
pri čemu ćemo \tilde{r} \geq 0.01m mijenjati, izrađen od čelika (\rho=7850{\text{k}g/m^{3}}, E=2.1\, 10^{11}Pa). Tipični položaj pri prve dvije kritične brzine prikazan je na Slikama 15, 16.


Slika 15: Položaj pri 1. kritičnoj brzini


Slika 16: Položaj pri 2. kritičnoj brzini


Zanimljivije je pogledati ovisnost (prve) kritične brzine o radijusu zadebljanog dijela \tilde{r}. Vrijednosti su dane u Tablici 3.


\tilde{r}[m] 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
\gamma_{\scriptsize{\mbox{crit}}}[rpm] 5533 4466 3294 2560 2082 1751 1509 1325 1180 1064
Tablica 3: Kritične brzine u okretajima u minuti u ovisnosti o radijusu \tilde{r} proširenog dijela štapa

Vidimo da kritična brzina opada s povećanjem radijusa, tj. sustav postaje nestabilniji, kao što je i očekivano.

Primjer 5.[Svrdlo, povećavamo kontaktnu silu] U ovom primjeru analiziramo promjenu kritične brzine pri promjeni kontaktne sile. Pri tome imamo u vidu primjer bušilice. Svrdlo na strani bušilice smatramo ukliještenim, dok na vrhu svrdla postavljamo rubni uvjet šarnira
(17)
u(0)= u''(0)= u(L)= u'(L)=0.
Duljina svrdla u ovom primjeru je L=0.6m, promjer je 2R=0.01m, dok je svrdlo napravljeno od čelika. I opet je zanimljivo pogledati kako se mijenja kritična brzina, ovaj put u ovisnosti o kontaktnoj sili \sigma, Tablica 4. Negativan predznak znači da je sila kompresija.


\sigma[N] 0 -1000 -2000 -3000 -4000 -5000 -5600 -5700 -5780 -5781
\gamma_{\scriptsize{\mbox{crit}}}[rpm] 5296 4826 4301 3697 2966 1970 951 638 86 49
Tablica 4: Kritične brzine u okretajima u minuti u ovisnosti o kompresivnoj kontaktnoj sili \sigma

Iz rezultata vidimo da kako povećavamo iznos sile, kritična brzina opada, a rezultati sugeriraju da kritična brzina ide u 0, kako i teorija predviđa, i dovodi do Eulerove nestabilnosti ili bucklinga, tj. nestabilnosti i bez vrtnje. Položaj pri prve dvije kritične brzine dan je na slikama 17 i 18.


Slika 17: Položaj pri 1. kritičnoj brzini


Slika 18: Položaj pri 2. kritičnoj brzini


Primjer 6.[Svrdlo, mijenjamo duljinu] U ovom primjeru ponovo promatramo svrdlo, pa su nam rubni uvjeti isti kao u (17). Stoga su položaji pri kritičnim brzinama kao u prethodnom primjeru (slikama 17 i 18). Razlika je u tome da sada analiziramo ovisnost o duljini svrdla. Ovo je posebno važno npr. u geologiji, posebno kod naftnih bušotina. Promjer je i dalje 2R=0.01m, a svrdlo je čelično. Duljinu uzimamo od 1m do 40m. Kritične brzine su u Tablici 5.




L[m] 1 5 10 15 20 40
\gamma_{\scriptsize{\mbox{crit}}}[rpm] 1906 76 19 8 5 1
Tablica 5: Kritične brzine u okretajima u minuti u ovisnosti o duljini svrdla

Naravno, svrdlo promjera 1cm i duljine 40m nikad nismo vidjeli, a razlog se vidi u ovim rezultatima.

Primjer 7.[Osovina generatora] Jedan mogući izgled osovine generatora dan je na Slici 19.


Slika 19: Oblik osovine generatora

Duljina osovine je L=2.165m, presjeci su kružni maksimalnog promjera 2R=0.116m, a za materijal izrade uzet je čelik. Na slikama 20 i 21 su položaji osovine pri prve dvije kritične brzine.


Slika 20: Položaj pri 1. kritičnoj brzini


Slika 21: Položaj pri 2. kritičnoj brzini


Dobivene kritične brzine nalaze se u Tablici 6.




kritična brzina 1 2 3 4
\gamma_{\scriptsize{\mbox{i}}}[rpm] 492 2225 5654 62075
Tablica 6: kritične brzine u okretajima u minutu u ovisnosti o duljini svrdla

Ovi rezultati ipak nisu sasvim realni jer se provodi dodatna stabilizacija osovine koja nije uzeta u obzir u promatranom modelu.

7Eksperiment

Da bismo neki model upotrijebili u praksi, najčešće je nužna verifikacija modela s pomoću eksperimenta. Postavljanje i vođenje eksperimenta iznimno je zahtjevan zadatak te traži puno inženjerskog znanja i iskustva. Fizikalni model izrađen je i ispitivan na Katedri za turbostrojeve Fakulteta strojarstva i brodogradnje Sveučilišta u Zagrebu. Eksperiment su organizirali i vodili prof. dr. sc. Branimir Matijašević i prof. dr. sc. Zvonimir Guzović.


Slika 22: Postavka eksperimenta

Osnovni podaci eksperimenta su sljedeći:
\bullet duljina štapa: L=0.835m
\bullet promjer štapa: 2R=0.01m
\bullet uzdužna napetost: \sigma = 0N
\bullet materijal štapa: čelik
\bullet Youngov modul elastičnosti: E=2.1 \, 10^{11}Pa
\bullet gustoća: \varrho = 7800 \frac{\text{kg}}{\text{m}^{3}}
Obavljena su dva eksperimenta — jedan za ukliještene krajeve i jedan za krajeve učvršćene šarnirom. Za određivanje kritične brzine korišten je stroboskop (Slika 23),


Slika 23: Stroboskop

dok je za rotaciju štapa korišten elektromotor (Slika 24).


Slika 24: Elektromotor

Ograničenja korištenih uređaja omogućila su određivanje samo prvih dviju kritičnih brzina.

7.11. eksperiment: šarnir

U slučaju ovih rubnih uvjeta, matematički model te eksperiment daju rezultate u Tablici 7.




  \gamma_{1}[rpm] \gamma_{1}[rpm] \gamma_{1}[rpm] \gamma_{1}[rpm]
Model 1754 7015 15788 28076
Eksperiment 1744 6830 - -
Tablica 7: Rezultati dobiveni iz matematičkog modela i eksperimenta

Iz rezultata vidimo da je razlika za prvu kritičnu brzinu ispod 1%, a za drugu ispod 3%.

7.22. eksperiment: ukliješteni krajevi

U slučaju ovih rubnih uvjeta, matematički model te eksperiment daju rezultate u Tablici 8.




  \gamma_{1}[rpm] \gamma_{1}[rpm] \gamma_{1}[rpm] \gamma_{1}[rpm]
Model 3975 10959 21490 35536
Eksperiment 1855 3480 - -
Tablica 8: Rezultati dobiveni iz matematičkog modela i eksperimenta

Rezultati eksperimenta jako odudaraju od onih dobivenih iz modela. Naravno, u ovom slučaju delikatna je implementacija rubnih uvjeta. Za primijetiti je da je 1. kritična brzina zapravo vrlo bliska onoj iz 1. eksperimenta.

Bibliografija
[1] Aganović, I.; Tambača, J.: On the stability of rotating rods and plates, ZAMM 81 (2001), 733. – 742.
[2] Aganović, I.; Veselić, K.: Jednadžbe matematičke fizike, školska knjiga, Zagreb, 1985.
[3] Kurepa, S.: Funkcionalna analiza, školska knjiga, Zagreb, 1994.
[4] Landau, L.; Lifchitz, E.: Theory of Elasticity. Pergamon Press, 1970.
[5] Tokmačić, J.: Stabilnost rotirajućeg štapa, Diplomski rad, Sveučilište u Zagrebu, PMF-Matematički odjel, Zagreb, prosinac 2005.
[6] Trabucho, L.; Viaño, J. M.: Mathematical modelling of rods. In: Ciarlet, P. G.; Lions, J. L. (eds.): Handbook of Numerical Analysis, Vol. 4. North–Holland, 1996.
[7] Tutek, Z.; Aganović, I.: A justification of the one-dimensional linear model of elastic beam, Math. Models Appl. Sci. 8 (1986), 502. – 515.
[8] Tutek, Z.; Vrdoljak, M: Obične diferencijalne jednadžbe, skripta PMF-Matematički odjel, 2009.