Модальний аналіз¶

Узагальнена задача на власні значення¶

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

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

Тут $u$ — узагальнений вектор переміщень, $M$ — матриця мас, а $K$ — матриця жорсткості. Нехай $\omega$ — власна кутова частота, $a$ і $b$ — довільні сталі, що не дорівнюють нулю одночасно, а $x$ — вектор. Визначимо функцію

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

Її друга похідна дорівнює

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

Підставивши рівняння $\eqref{eq:2.3.2}$ і рівняння $\eqref{eq:2.3.3}$ у рівняння $\eqref{eq:2.3.1}$, отримаємо

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

Отже, отримуємо наведене вище рівняння.

Для нетривіального коливання вираз $a \sin \omega t + b \cos \omega t$ не є тотожно нульовим.

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

Тому виконується наведене вище співвідношення. Поклавши $\lambda=\omega^2$, отримаємо

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

Отже, отримуємо наведене вище рівняння.

Скаляр $\lambda$ називають власним значенням, вектор $x$ — власним вектором, а задачу, виражену рівнянням $\eqref{eq:2.3.5}$, — узагальненою задачею на власні значення.

Власну кутову частоту $\omega$ визначають із власного значення $\lambda=\omega^2$, а відповідний власний вектор $x$ задає форму коливань.

Приклад системи з багатьма ступенями свободи для вільних коливань без демпфування

Рисунок 2.3.1 Приклад системи з багатьма ступенями свободи для вільних коливань без демпфування

Властивості матриць і припущення¶

Для узагальненої задачі на власні значення $Kx=\lambda Mx$, отриманої в попередньому розділі, у цьому посібнику припускається симетричність матриць. Для комплексних матриць це відповідає ермітовим матрицям, а для дійсних — симетричним матрицям.

Нехай $k_{ij}$ — елемент $ij$ матриці $K$. Ермітовість виражається співвідношенням

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

у наведеному вище рівнянні, де $\bar{k}_{ji}$ — комплексно спряжене до $k_{ji}$. Для дійсних матриць це співвідношення набуває вигляду $k_{ij}=k_{ji}$.

Дійсна симетрична матриця $H$ є додатно визначеною, якщо для кожного ненульового вектора $x$

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

виконується наведена вище нерівність. У цьому разі всі власні значення матриці $H$ додатні.

У структурній задачі на власні значення матрицю мас $M$ зазвичай вважають додатно визначеною. Матриця жорсткості $K$, однак, залежно від умов закріплення може бути додатно напіввизначеною та мати нульові власні значення, що відповідають модам абсолютно твердого тіла.

Обернена ітерація зі зсувом¶

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

Для зсуву $\sigma$, якщо $-\sigma$ не є власним значенням, а $K+\sigma M$ є невиродженою, рівняння $\eqref{eq:2.3.5}$ можна перетворити так.

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

Це перетворення не змінює власний вектор $x$ і відображає власне значення $\lambda$ у $1/(\lambda+\sigma)$.

Отже, що ближче $\lambda$ до $-\sigma$, то більшим є модуль перетвореного власного значення. У структурній задачі на власні значення $\lambda \geq 0$ і $\sigma \geq 0$, тому найнижче власне значення відображається у власне значення з найбільшим модулем. Використовуючи цю властивість і застосовуючи до рівняння $\eqref{eq:2.3.8}$ ітераційний метод, який спочатку збігається до власних значень із великим модулем, можна ефективно отримувати власні значення нижчих порядків у зростаючому порядку.

Цей метод називають оберненою ітерацією зі зсувом.

Для аналізів із закріпленнями FrontISTR задає $\sigma = 0$, і рівняння $\eqref{eq:2.3.8}$ набуває вигляду $K^{-1} M x = \frac{1}{\lambda} x$, тобто стає оберненою ітерацією без зсуву. У незакріпленому аналізі вільного тіла матриця $K$ є виродженою, оскільки має нульові власні значення, що відповідають модам абсолютно твердого тіла; додатне значення $\sigma$ регуляризує $K+\sigma M$. Значення $\sigma$ задається параметром SIGMA в !EIGEN.

Метод Ланцоша¶

Обґрунтування вибору (порівняння з методом Якобі)¶

Метод Якобі є добре відомим класичним методом.

Він ефективний, коли матриця мала й щільна. Однак, оскільки матриці, які обробляє HEC-MW, великі й розріджені, метод Якобі не використовується; натомість застосовується ітераційний метод Ланцоша.

Алгоритм і особливості¶

Цей метод, запропонований К. Ланцошем у 1950-х роках, є алгоритмом тридіагоналізації матриці та має такі особливості.

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

Метод Ланцоша починається з початкового вектора, послідовно будує ортогональні вектори та отримує базис підпростору Крилова.

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

Геометричний зміст (підпростір Крилова)¶

Застосувавши до рівняння $\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}\]

задачу можна переписати у вигляді

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

Отже, отримуємо наведене вище рівняння.

Для довільного ненульового вектора $q_0$

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

простір, натягнутий на ці вектори,

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

називають підпростором Крилова.

Метод Ланцоша послідовно будує базис цього підпростору Крилова.

FrontISTR використовує такий скалярний добуток, пов’язаний із матрицею мас $M$,

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

для ортонормування базису. Скалярний добуток $\langle x,y\rangle$ і норма $\|x\|$, показані на наступних рисунках, в обчисленнях FrontISTR інтерпретуються відповідно як цей $M$-скалярний добуток і відповідна $M$-норма, визначена співвідношенням:

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

Застосуйте лінійне перетворення, визначене матрицею $A$, до довільного вектора $q_0$, як показано на рисунку 2.3.2.

$Лінійне перетворення $q_0$ матрицею $A$$

Рисунок 2.3.2 Лінійне перетворення $q_0$ матрицею $A$

Перетворений вектор ортогоналізується в просторі, натягнутому на нього та початковий вектор. Зокрема, виконується ортогоналізація Грама—Шмідта, як показано на рисунку 2.3.3. Позначимо отриманий вектор через $r_1$; його нормування дає $q_1$.

$Вектор $q_1$, ортогональний до $q_0$$

Рисунок 2.3.3 Вектор $q_1$, ортогональний до $q_0$

Тією самою процедурою отримуємо $q_2$ з $q_1$. При цьому $q_2$ ортогональний і до $q_1$, і до $q_0$, як показано на рисунку 2.3.4.

$Вектор $q_2$, ортогональний до $q_1$ і $q_0$$

Рисунок 2.3.4 Вектор $q_2$, ортогональний до $q_1$ і $q_0$

Таким чином, метод Ланцоша послідовно будує ортонормований базис підпростору Крилова. Теоретично симетрія задачі на власні значення дає змогу подати цю ітерацію як тричленне рекурентне співвідношення з використанням останніх векторів базису.

Однак реалізація FrontISTR повторно ортогоналізує кожен вектор відносно вже отриманих векторів базису Ланцоша за допомогою $M$-скалярного добутку, щоб зменшити втрату ортогональності через арифметику зі скінченною точністю.

Тридіагоналізація¶

В ітерації Ланцоша FrontISTR вектори базису ортонормуються відносно $M$-скалярного добутку, описаного в попередньому розділі, тому

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

виконується наведене вище співвідношення.

З використанням матриці $A$ з рівняння $\eqref{eq:2.3.10}$ теоретична ітерація Ланцоша подається таким тричленним рекурентним співвідношенням:

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

Спочатку $\alpha_i$ визначається так.

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

Далі визначимо попередній вектор нев’язки як

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

як показано вище.

Щоб зменшити втрату ортогональності через арифметику зі скінченною точністю, реалізація FrontISTR повторно ортогоналізує $\tilde{r}_{i+1}$ відносно вже отриманих векторів базису Ланцоша за допомогою $M$-скалярного добутку. Якщо $r_{i+1}$ — вектор нев’язки після повторної ортогоналізації, то

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

виконуються наведені вище співвідношення.

Об’єднаємо $m$ векторів базису, отриманих ітерацією Ланцоша, у

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

Тоді після скінченної кількості ітерацій Ланцоша

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

виконується наведене вище співвідношення.

Тут $e_m$ — $m$-вимірний одиничний вектор, $m$-та компонента якого дорівнює 1, а всі інші компоненти дорівнюють нулю, і

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

є симетричною тридіагональною матрицею.

Отже, власні значення початкової великомасштабної задачі можна наближено визначити, обчисливши власні значення тридіагональної матриці $T_m$.

Пов’язані теми¶

Аналіз частотної характеристики — Частотна характеристика з використанням результатів модального аналізу
Типи аналізу — Функціональний огляд модального аналізу