Přeskočit obsah

Modální analýza

Zobecněná úloha vlastních čísel

Při analýze volného kmitání kontinua se konstrukce prostorově diskretizuje a modeluje jako soustava s více stupni volnosti a soustředěnými hmotami, jak je znázorněno na obrázku 2.3.1. Pro úlohu netlumeného volného kmitání má řídicí rovnice (pohybová rovnice) následující tvar.

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

Zde je \(u\) zobecněný vektor posunutí, \(M\) matice hmotnosti a \(K\) matice tuhosti. Nechť \(\omega\) je vlastní úhlová frekvence, \(a\) a \(b\) jsou libovolné konstanty, které nejsou současně nulové, a \(x\) je vektor. Definujme funkci

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

Její druhá derivace je

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

Dosazením rovnice \(\eqref{eq:2.3.2}\) a rovnice \(\eqref{eq:2.3.3}\) do rovnice \(\eqref{eq:2.3.1}\) dostaneme

\[\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}\]

Tím dostaneme výše uvedenou rovnici.

Pro netriviální kmitání není \(a \sin \omega t + b \cos \omega t\) identicky rovno nule.

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

Proto platí výše uvedený vztah. Položíme-li \(\lambda=\omega^2\), dostaneme

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

Tím dostaneme výše uvedenou rovnici.

Skalár \(\lambda\) se nazývá vlastní číslo, vektor \(x\) vlastní vektor a úloha vyjádřená rovnicí \(\eqref{eq:2.3.5}\) zobecněná úloha vlastních čísel.

Vlastní úhlová frekvence \(\omega\) se získá z vlastního čísla \(\lambda=\omega^2\) a odpovídající vlastní vektor \(x\) představuje vlastní tvar.

Příklad soustavy s více stupni volnosti pro netlumené volné kmitání

Obrázek 2.3.1 Příklad soustavy s více stupni volnosti pro netlumené volné kmitání

Vlastnosti matic a předpoklady

Pro zobecněnou úlohu vlastních čísel \(Kx=\lambda Mx\) získanou v předchozím oddílu tato příručka předpokládá, že matice jsou symetrické. U komplexních matic tomu odpovídají hermitovské matice, u reálných matic symetrické matice.

Nechť \(k_{ij}\) je prvek \(ij\) matice \(K\). Hermitovská vlastnost je vyjádřena vztahem

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

uvedeným výše, kde \(\bar{k}_{ji}\) je komplexně sdružená hodnota k \(k_{ji}\). Pro reálné matice se tento vztah změní na \(k_{ij}=k_{ji}\).

Reálná symetrická matice \(H\) je pozitivně definitní, jestliže pro každý nenulový vektor \(x\)

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

platí výše uvedená nerovnost. V takovém případě jsou všechna vlastní čísla matice \(H\) kladná.

Ve strukturní úloze vlastních čísel se matice hmotnosti \(M\) obvykle považuje za pozitivně definitní. Matice tuhosti \(K\) však může být v závislosti na okrajových podmínkách pozitivně semidefinitní a může mít nulová vlastní čísla odpovídající módům tuhého tělesa.

Inverzní iterace s posunem

V praktické strukturní analýze metodou konečných prvků zpravidla nejsou zapotřebí všechna vlastní čísla; často postačí nanejvýš několik vlastních čísel nízkého řádu. HEC-MW je určen pro rozsáhlé úlohy, v nichž jsou matice velké a velmi řídké, s mnoha nulovými prvky. S ohledem na tyto vlastnosti je proto důležité efektivně počítat vlastní čísla módů nízkého řádu.

Pro posun \(\sigma\), pokud \(-\sigma\) není vlastním číslem a \(K+\sigma M\) je regulární, lze rovnici \(\eqref{eq:2.3.5}\) převést následovně.

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

Tato transformace ponechá vlastní vektor \(x\) beze změny a zobrazí vlastní číslo \(\lambda\) na \(1/(\lambda+\sigma)\).

Čím blíže je tedy \(\lambda\) hodnotě \(-\sigma\), tím větší je absolutní hodnota transformovaného vlastního čísla. Ve strukturní úloze vlastních čísel platí \(\lambda \geq 0\) a \(\sigma \geq 0\), takže nejnižší vlastní číslo se zobrazí na vlastní číslo s největší absolutní hodnotou. Využitím této vlastnosti a použitím iterační metody, která má tendenci nejprve konvergovat k vlastním číslům s velkou absolutní hodnotou, na rovnici \(\eqref{eq:2.3.8}\) lze efektivně získat vlastní čísla nízkého řádu ve vzestupném pořadí.

Tato metoda se nazývá inverzní iterace s posunem.

U analýz s vazbami nastaví FrontISTR \(\sigma = 0\) a rovnice \(\eqref{eq:2.3.8}\) přejde na \(K^{-1} M x = \frac{1}{\lambda} x\), tedy na inverzní iteraci bez posunu. Při nevázané analýze volné konstrukce je \(K\) singulární, protože má nulová vlastní čísla odpovídající módům tuhého tělesa; přiřazením kladné hodnoty \(\sigma\) se matice \(K+\sigma M\) regularizuje. Hodnota \(\sigma\) se zadává parametrem SIGMA v !EIGEN.

Lanczosova metoda

Důvod použití (porovnání s Jacobiho metodou)

Jacobiho metoda je dobře známá klasická metoda.

Je účinná, když je matice malá a hustá. Protože jsou však matice zpracovávané systémem HEC-MW velké a řídké, Jacobiho metoda se nepoužívá; namísto ní se používá iterační Lanczosova metoda.

Algoritmus a vlastnosti

Tato metoda, navržená C. Lanczosem v 50. letech 20. století, je algoritmem pro tridiagonalizaci matice a má následující vlastnosti.

  • Je to iterační metoda a umožňuje provádět výpočet při zachování řídkosti matice.
  • Její operace tvoří především součiny matice a vektoru, takže je vhodná pro paralelizaci.
  • Je vhodná pro geometrickou dekompozici oblasti související se sítěmi konečných prvků.
  • Efektivní výpočet lze provést omezením počtu hledaných vlastních čísel nebo rozsahu módů.

Lanczosova metoda vychází z počátečního vektoru, postupně vytváří ortogonální vektory a získává bázi Krylovova podprostoru.

Při aritmetice s konečnou přesností mohou zaokrouhlovací chyby způsobit ztrátu ortogonality vektorů. Aby se tento vliv potlačil, implementace FrontISTR provádí opětovnou ortogonalizaci každého vektoru vůči již získaným vektorům Lanczosovy báze.

Geometrický význam (Krylovův podprostor)

Použitím následujících transformací proměnných na rovnici \(\eqref{eq:2.3.8}\)

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

lze úlohu přepsat do tvaru

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

Tím dostaneme výše uvedenou rovnici.

Pro libovolný nenulový vektor \(q_0\)

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

se prostor generovaný těmito vektory

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

nazývá Krylovův podprostor.

Lanczosova metoda postupně konstruuje bázi tohoto Krylovova podprostoru.

FrontISTR používá následující skalární součin spojený s maticí hmotnosti \(M\)

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

k ortonormalizaci báze. Skalární součin \(\langle x,y\rangle\) a norma \(\|x\|\) znázorněné na následujících obrázcích se ve výpočtech FrontISTR interpretují jako tento \(M\)-skalární součin a odpovídající \(M\)-norma

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

v uvedeném pořadí.

Použijte lineární transformaci definovanou maticí \(A\) na libovolný vektor \(q_0\), jak je znázorněno na obrázku 2.3.2.

Lineární transformace vektoru \(q_0\) maticí \(A\)

Obrázek 2.3.2 Lineární transformace vektoru \(q_0\) maticí \(A\)

Transformovaný vektor se ortogonalizuje v prostoru generovaném jím a původním vektorem. Konkrétně se provede Gramova-Schmidtova ortogonalizace, jak je znázorněno na obrázku 2.3.3. Označme výsledný vektor \(r_1\); jeho normalizací získáme \(q_1\).

Vektor \(q_1\) kolmý k \(q_0\)

Obrázek 2.3.3 Vektor \(q_1\) kolmý k \(q_0\)

Stejným postupem získáme \(q_2\) z \(q_1\). Vektor \(q_2\) je přitom kolmý k \(q_1\) i \(q_0\), jak je znázorněno na obrázku 2.3.4.

Vektor \(q_2\) kolmý k \(q_1\) a \(q_0\)

Obrázek 2.3.4 Vektor \(q_2\) kolmý k \(q_1\) a \(q_0\)

Lanczosova metoda tak postupně konstruuje ortonormální bázi Krylovova podprostoru. Teoreticky symetrie úlohy vlastních čísel umožňuje vyjádřit tuto iteraci jako tříčlennou rekurenci využívající nejnovější vektory báze.

Implementace FrontISTR však provádí opětovnou ortogonalizaci každého vektoru vůči již získaným vektorům Lanczosovy báze pomocí \(M\)-skalárního součinu, aby potlačila ztrátu ortogonality způsobenou aritmetikou s konečnou přesností.

Tridiagonalizace

V Lanczosově iteraci systému FrontISTR se vektory báze ortonormalizují vzhledem k \(M\)-skalárnímu součinu popsanému v předchozím oddílu, takže

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

platí výše uvedený vztah.

Při použití matice \(A\) z rovnice \(\eqref{eq:2.3.10}\) má teoretická Lanczosova iterace tvar

\[\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}\]

vyjádřený výše uvedenou tříčlennou rekurencí.

Nejprve se \(\alpha_i\) definuje následovně.

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

Poté se předběžné reziduum definuje jako

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

jak je uvedeno výše.

Aby se potlačila ztráta ortogonality způsobená aritmetikou s konečnou přesností, implementace FrontISTR znovu ortogonalizuje \(\tilde{r}_{i+1}\) vůči již získaným vektorům Lanczosovy báze pomocí \(M\)-skalárního součinu. Je-li \(r_{i+1}\) reziduum po opětovné ortogonalizaci, pak

\[\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}\]

platí výše uvedené vztahy.

Sestavme \(m\) vektorů báze získaných Lanczosovou iterací do matice

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

Potom po konečném počtu Lanczosových iterací

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

platí výše uvedený vztah.

Zde je \(e_m\) \(m\)-rozměrný jednotkový vektor, jehož \(m\)-tá složka je 1 a ostatní složky jsou nulové, a

\[\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}\]

je symetrická tridiagonální matice.

Vlastní čísla původní rozsáhlé úlohy vlastních čísel lze tedy aproximovat výpočtem vlastních čísel tridiagonální matice \(T_m\).

Související témata