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

Обобщённая задача на собственные значения¶

При анализе свободных колебаний сплошной среды выполняется пространственная дискретизация, после чего она моделируется как система со многими степенями свободы и сосредоточенными массами, показанная на рис. 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$ в настоящем руководстве предполагается симметричность рассматриваемых матриц. Для комплексных матриц это соответствует эрмитовым матрицам, а для вещественных — симметричным.

Обозначим элемент $ij$ матрицы $K$ через $k_{ij}$. Эрмитовость выражается соотношением

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

соответственно.

К произвольному вектору $q_0$ применяется линейное преобразование, заданное матрицей $A$, как показано на рис. 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$.

Связанные разделы¶

Анализ частотного отклика — частотный отклик с использованием результатов модального анализа
Типы анализа — обзор возможностей модального анализа