Analiza modalna¶

Uogólnione zagadnienie własne¶

W analizie drgań swobodnych ośrodka ciągłego konstrukcję dyskretyzuje się przestrzennie i modeluje jako układ o wielu stopniach swobody z masami skupionymi, jak pokazano na rysunku 2.3.1. Dla zagadnienia swobodnych drgań nietłumionych równanie rządzące (równanie ruchu) ma następującą postać.

\[\begin{equation} M \ddot{u} + K u = 0 \label{eq:2.3.1} \end{equation}\]

Tutaj $u$ jest uogólnionym wektorem przemieszczeń, $M$ jest macierzą masy, a $K$ jest macierzą sztywności. Niech $\omega$ oznacza naturalną częstość kołową, $a$ i $b$ będą dowolnymi stałymi, które nie są jednocześnie równe zeru, a $x$ będzie wektorem. Zdefiniujmy funkcję

\[\begin{equation} u(t) = (a \sin \omega t + b \cos \omega t ) x \label{eq:2.3.2} \end{equation}\]

Jej druga pochodna ma postać

\[\begin{equation} \ddot{u}(t) = -\omega^2 (a \sin \omega t + b \cos \omega t) x \label{eq:2.3.3} \end{equation}\]

Po podstawieniu równań $\eqref{eq:2.3.2}$ i $\eqref{eq:2.3.3}$ do równania $\eqref{eq:2.3.1}$ otrzymuje się następującą zależność:

\[\begin{equation} M \ddot{u} + K u = (a \sin \omega t + b \cos \omega t) (K-\omega^2 M) x = 0 \label{eq:2.3.4} \end{equation}\]

Powyższa zależność jest wynikiem tego podstawienia.

Dla drgania nietrywialnego wyrażenie $a \sin \omega t + b \cos \omega t$ nie jest tożsamościowo równe zeru, dlatego zachodzi zależność:

\[ (K-\omega^2M)x=0 \]

Po przyjęciu $\lambda=\omega^2$ zależność tę można zapisać w postaci:

\[\begin{equation} K x = \lambda M x \label{eq:2.3.5} \end{equation}\]

Otrzymuje się w ten sposób powyższe równanie.

Skalar $\lambda$ nazywa się wartością własną, wektor $x$ — wektorem własnym, a zagadnienie opisane równaniem $\eqref{eq:2.3.5}$ — uogólnionym zagadnieniem własnym.

Naturalną częstość kołową $\omega$ wyznacza się z wartości własnej $\lambda=\omega^2$, a odpowiadający jej wektor własny $x$ przedstawia postać drgań.

Przykład układu o wielu stopniach swobody dla nietłumionych drgań swobodnych

Rysunek 2.3.1 Przykład układu o wielu stopniach swobody dla nietłumionych drgań swobodnych

Własności i założenia dotyczące macierzy¶

Dla uogólnionego zagadnienia własnego $Kx=\lambda Mx$, otrzymanego w poprzednim rozdziale, w niniejszym podręczniku zakłada się symetrię macierzy. W przypadku macierzy zespolonych odpowiada to macierzom hermitowskim, a w przypadku macierzy rzeczywistych — macierzom symetrycznym.

Niech $k_{ij}$ oznacza element $ij$ macierzy $K$. Własność hermitowską wyraża następująca zależność:

\[\begin{equation} k_{ij} = \bar{k}_{ji} \label{eq:2.3.6} \end{equation}\]

W powyższej zależności $\bar{k}_{ji}$ oznacza sprzężenie zespolone elementu $k_{ji}$. Dla macierzy rzeczywistych zależność ta przyjmuje postać $k_{ij}=k_{ji}$.

Rzeczywista macierz symetryczna $H$ jest dodatnio określona, jeżeli dla każdego niezerowego wektora $x$ zachodzi nierówność:

\[\begin{equation} x^{t} H x > 0 \label{eq:2.3.7} \end{equation}\]

W takim przypadku wszystkie wartości własne macierzy $H$ są dodatnie.

W konstrukcyjnym zagadnieniu własnym macierz masy $M$ jest zwykle traktowana jako dodatnio określona. Macierz sztywności $K$ może natomiast, zależnie od warunków podparcia, być dodatnio półokreślona i mieć zerowe wartości własne odpowiadające modom bryły sztywnej.

Odwrotna iteracja z przesunięciem¶

W praktycznej analizie konstrukcji metodą elementów skończonych zwykle nie są potrzebne wszystkie wartości własne; często wystarcza najwyżej kilka wartości własnych najniższych rzędów. HEC-MW jest przeznaczony do zagadnień wielkoskalowych, w których macierze są duże i bardzo rzadkie, z wieloma elementami zerowymi. Dlatego ważne jest efektywne obliczanie wartości własnych modów niskiego rzędu z uwzględnieniem tych cech.

Dla przesunięcia $\sigma$, jeśli $-\sigma$ nie jest wartością własną, a $K+\sigma M$ jest nieosobliwa, równanie $\eqref{eq:2.3.5}$ można przekształcić następująco.

\[\begin{equation} (K + \sigma M)^{-1} M x = \frac{1}{\lambda+\sigma} x \label{eq:2.3.8} \end{equation}\]

Przekształcenie to nie zmienia wektora własnego $x$, a wartość własną $\lambda$ odwzorowuje na $1/(\lambda+\sigma)$.

W konsekwencji im bliżej $\lambda$ znajduje się wartości $-\sigma$, tym większy jest moduł przekształconej wartości własnej. W konstrukcyjnym zagadnieniu własnym $\lambda \geq 0$ i $\sigma \geq 0$, dlatego najniższa wartość własna zostaje odwzorowana na wartość własną o największym module. Wykorzystując tę własność i stosując do równania $\eqref{eq:2.3.8}$ metodę iteracyjną, która ma tendencję do wcześniejszej zbieżności do wartości własnych o dużym module, można efektywnie wyznaczać wartości własne niskich rzędów w kolejności rosnącej.

Metodę tę nazywa się odwrotną iteracją z przesunięciem.

W analizach z więzami FrontISTR przyjmuje $\sigma = 0$, a równanie $\eqref{eq:2.3.8}$ przyjmuje postać $K^{-1} M x = \frac{1}{\lambda} x$, czyli odwrotnej iteracji bez przesunięcia. W analizie swobodnej bez więzów macierz $K$ jest osobliwa, ponieważ ma zerowe wartości własne odpowiadające modom bryły sztywnej; nadanie $\sigma$ wartości dodatniej powoduje regularyzację macierzy $K+\sigma M$. Wartość $\sigma$ określa parametr SIGMA w !EIGEN.

Metoda Lanczosa¶

Uzasadnienie wyboru (porównanie z metodą Jacobiego)¶

Metoda Jacobiego jest dobrze znaną metodą klasyczną.

Jest skuteczna, gdy macierz jest mała i gęsta. Ponieważ jednak macierze obsługiwane przez HEC-MW są duże i rzadkie, metoda Jacobiego nie jest stosowana; zamiast niej używa się iteracyjnej metody Lanczosa.

Algorytm i cechy¶

Metoda ta, zaproponowana przez C. Lanczosa w latach 50. XX wieku, jest algorytmem trójdiagonalizacji macierzy i ma następujące cechy.

Jest metodą iteracyjną i umożliwia prowadzenie obliczeń przy zachowaniu rzadkiej postaci macierzy.
Jej operacje obejmują głównie mnożenie macierzy przez wektor, dzięki czemu dobrze nadaje się do zrównoleglania.
Dobrze nadaje się do geometrycznego podziału domeny związanego z siatkami elementów skończonych.
Umożliwia efektywne obliczenia przez ograniczenie liczby wyznaczanych wartości własnych lub zakresu modów.

Metoda Lanczosa rozpoczyna się od wektora początkowego, kolejno tworzy wektory ortogonalne i wyznacza bazę podprzestrzeni Kryłowa.

W arytmetyce o skończonej precyzji błędy zaokrągleń mogą powodować utratę ortogonalności wektorów. Aby ograniczyć ten efekt, implementacja FrontISTR ponownie ortogonalizuje każdy wektor względem już wyznaczonych wektorów bazy Lanczosa.

Interpretacja geometryczna (podprzestrzeń Kryłowa)¶

Wprowadźmy do równania $\eqref{eq:2.3.8}$ następujące przekształcenia zmiennych:

\[ A = (K + \sigma M)^{-1} M \]

\[\begin{equation} \frac{1}{\lambda+\sigma}= \zeta \label{eq:2.3.9} \end{equation}\]

Po ich zastosowaniu zagadnienie można zapisać w postaci:

\[\begin{equation} A x = \zeta x \label{eq:2.3.10} \end{equation}\]

W ten sposób zagadnienie przyjmuje przedstawioną wyżej postać.

Dla dowolnego niezerowego wektora $q_0$

\[ q_0,\ Aq_0,\ A^2q_0,\ldots,A^{m-1}q_0 \]

przestrzeń rozpięta przez te wektory

\[ \mathcal{K}_m(A,q_0) = \operatorname{span} \{q_0,Aq_0,A^2q_0,\ldots,A^{m-1}q_0\} \]

nazywa się podprzestrzenią Kryłowa.

Metoda Lanczosa kolejno konstruuje bazę tej podprzestrzeni Kryłowa.

Do ortonormalizacji bazy FrontISTR wykorzystuje następujący iloczyn skalarny związany z macierzą masy $M$:

\[ \langle x,y\rangle_M = x^T M y \]

W obliczeniach FrontISTR iloczyn skalarny $\langle x,y\rangle$ i norma $\|x\|$, pokazane na kolejnych rysunkach, są interpretowane odpowiednio jako powyższy iloczyn skalarny względem $M$ oraz odpowiadająca mu norma względem $M$, zdefiniowana wzorem:

\[ \|x\|_M=\sqrt{x^T M x} \]

W dalszych rysunkach oba symbole mają właśnie takie znaczenie.

Zastosuj przekształcenie liniowe określone przez macierz $A$ do dowolnego wektora $q_0$, jak pokazano na rysunku 2.3.2.

$Przekształcenie liniowe wektora $q_0$ przez macierz $A$$

Rysunek 2.3.2 Przekształcenie liniowe wektora $q_0$ przez macierz $A$

Przekształcony wektor jest ortogonalizowany w przestrzeni rozpiętej przez niego i wektor pierwotny. W szczególności wykonuje się ortogonalizację Grama-Schmidta, jak pokazano na rysunku 2.3.3. Niech otrzymany wektor będzie oznaczony przez $r_1$; po jego normalizacji otrzymuje się $q_1$.

$Wektor $q_1$ ortogonalny do $q_0$$

Rysunek 2.3.3 Wektor $q_1$ ortogonalny do $q_0$

Tę samą procedurę stosuje się do wyznaczenia $q_2$ z $q_1$. Wektor $q_2$ jest wówczas ortogonalny zarówno do $q_1$, jak i do $q_0$, jak pokazano na rysunku 2.3.4.

$Wektor $q_2$ ortogonalny do $q_1$ i $q_0$$

Rysunek 2.3.4 Wektor $q_2$ ortogonalny do $q_1$ i $q_0$

W ten sposób metoda Lanczosa kolejno konstruuje ortonormalną bazę podprzestrzeni Kryłowa. Teoretycznie symetria zagadnienia własnego pozwala zapisać tę iterację jako trójczłonową relację rekurencyjną wykorzystującą najnowsze wektory bazy.

Implementacja FrontISTR ponownie ortogonalizuje jednak każdy wektor względem wcześniej otrzymanych wektorów bazy Lanczosa przy użyciu iloczynu skalarnego względem $M$, aby ograniczyć utratę ortogonalności spowodowaną arytmetyką o skończonej precyzji.

Trójdiagonalizacja¶

W iteracji Lanczosa programu FrontISTR wektory bazy są ortonormalizowane w sensie opisanego w poprzednim podrozdziale iloczynu skalarnego względem $M$. Zachodzi zatem relacja:

\[ q_i^T M q_j = \delta_{ij} \]

Relacja ta wyraża ortonormalność bazy.

Przy użyciu macierzy $A$ z równania $\eqref{eq:2.3.10}$ teoretyczna iteracja Lanczosa ma postać:

\[\begin{equation} A q_i = \beta_i q_{i-1} + \alpha_i q_i + \beta_{i+1} q_{i+1} \label{eq:2.3.11} \end{equation}\]

Powyższe równanie jest trójczłonową relacją rekurencyjną.

Najpierw współczynnik $\alpha_i$ definiuje się wzorem:

\[ \alpha_i = q_i^T M A q_i \]

Następnie tymczasowy wektor resztowy definiuje się wzorem:

\[ \tilde{r}_{i+1} = Aq_i - \beta_i q_{i-1} - \alpha_i q_i \]

W ten sposób zostaje zdefiniowany tymczasowy wektor resztowy.

Aby ograniczyć utratę ortogonalności spowodowaną arytmetyką o skończonej precyzji, implementacja FrontISTR ponownie ortogonalizuje $\tilde{r}_{i+1}$ względem już uzyskanych wektorów bazy Lanczosa przy użyciu iloczynu skalarnego względem $M$. Jeżeli $r_{i+1}$ oznacza wektor resztowy po ponownej ortogonalizacji, otrzymuje się następujące zależności:

\[\begin{equation} \beta_{i+1} = \sqrt{r_{i+1}^T M r_{i+1}}, \qquad q_{i+1} = \frac{r_{i+1}}{\beta_{i+1}} \label{eq:2.3.12} \end{equation}\]

Obie wielkości są zatem wyznaczane przez powyższe wzory.

Niech poniższa macierz będzie utworzona z $m$ wektorów bazy otrzymanych w iteracji Lanczosa:

\[ Q_m=[q_0,q_1,\ldots,q_{m-1}] \]

Po skończonej liczbie iteracji Lanczosa zachodzi wówczas zależność:

\[\begin{equation} A Q_m = Q_m T_m + \beta_m q_m e_m^T \label{eq:2.3.13} \end{equation}\]

Otrzymuje się zatem powyższą relację macierzową.

Tutaj $e_m$ jest $m$-wymiarowym wektorem jednostkowym, którego $m$-ta składowa jest równa 1, a pozostałe składowe są równe zeru. Poniższa macierz ma postać:

\[\begin{equation} T_m= \begin{pmatrix} \alpha_0 & \beta_1 & & &\\ \beta_1 & \alpha_1 & \beta_2 & &\\ & \ddots & \ddots & \ddots &\\ & & \beta_{m-2} & \alpha_{m-2}& \beta_{m-1}\\ & & & \beta_{m-1} & \alpha_{m-1} \end{pmatrix} \label{eq:2.3.14} \end{equation}\]

Jest to symetryczna macierz trójdiagonalna.

Wartości własne pierwotnego wielkoskalowego zagadnienia własnego można zatem przybliżyć przez obliczenie wartości własnych macierzy trójdiagonalnej $T_m$.

Tematy powiązane¶

Analiza odpowiedzi częstotliwościowej — Odpowiedź częstotliwościowa z wykorzystaniem wyników analizy modalnej
Typy analiz — Przegląd funkcjonalny analizy modalnej