콘텐츠로 이동

모드 해석

일반화 고유값 문제

연속체의 자유 진동을 해석할 때는 공간 이산화를 수행하고, 그림 2.3.1과 같이 집중 질량을 갖는 다자유도계로 모델링한다. 감쇠가 없는 자유 진동 문제의 지배 방정식(운동 방정식)은 다음과 같다.

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

여기서 \(u\)는 일반화 변위 벡터, \(M\)은 질량 행렬, \(K\)는 강성 행렬이다. 고유 각진동수를 \(\omega\), 둘 다 0은 아닌 임의의 상수를 \(a\)와 \(b\), 벡터를 \(x\)라 하고 다음 함수를 정의한다.

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

이 함수의 2계 미분은 다음과 같다.

\[\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\)가 항등적으로 0이 아니므로

\[ (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}\)로 표현되는 문제를 일반화 고유값 문제라 한다.

고유값 \(\lambda=\omega^2\)에서 고유 각진동수 \(\omega\)를 구할 수 있으며, 대응하는 고유벡터 \(x\)는 진동 모드를 나타낸다.

감쇠가 없는 자유 진동 다자유도계의 예

그림 2.3.1　감쇠가 없는 자유 진동 다자유도계의 예

행렬의 성질과 가정

앞 절에서 얻은 일반화 고유값 문제 \(Kx=\lambda Mx\)에 대해 이 매뉴얼에서는 대상 행렬이 대칭이라고 가정한다. 복소 행렬에서는 에르미트 행렬에, 실수 행렬에서는 대칭 행렬에 해당한다.

행렬 \(K\)의 \(ij\) 성분을 \(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\)가 양의 정부호라는 것은 임의의 0이 아닌 벡터 \(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\)는 !EIGEN의 SIGMA로 지정한다.

Lanczos 법

채택 이유(Jacobi 법과의 비교)

고전적인 방법으로는 Jacobi 법이 잘 알려져 있다.

이 방법은 행렬의 크기가 작고 밀행렬일 때 효과적이다. 그러나 HEC-MW에서 다루는 행렬은 대규모 희소 행렬이므로 Jacobi 법을 사용하지 않고 Lanczos 반복법을 사용한다.

알고리즘과 특징

1950년대에 C. Lanczos가 제안한 이 방법은 행렬을 삼중대각화하는 알고리즘이며 다음과 같은 특징이 있다.

  • 반복법이므로 행렬의 희소성을 유지한 채 계산할 수 있다.
  • 연산의 대부분이 행렬-벡터 곱으로 구성되어 병렬화에 적합하다.
  • 유한요소 메시와 연계된 기하학적 영역 분할법에 적합하다.
  • 구할 고유값의 개수나 모드 범위를 제한하여 효율적으로 계산할 수 있다.

Lanczos 법은 초기 벡터에서 시작하여 직교 벡터를 차례로 구성하고 Krylov 부분공간의 기저를 구한다.

유한 정밀도 연산에서는 반올림 오차로 인해 벡터의 직교성이 손실될 수 있다. 이 영향을 억제하기 위해 FrontISTR 구현에서는 이미 구한 Lanczos 기저 벡터에 대해 재직교화를 수행한다.

기하학적 의미(Krylov 부분공간)

식 \(\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}\]

따라서 위 식을 얻는다.

임의의 0이 아닌 벡터 \(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\} \]

을 Krylov 부분공간이라 한다.

Lanczos 법은 이 Krylov 부분공간의 기저를 차례로 구성한다.

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

행렬 \(A\)에 의한 \(q_0\)의 선형 변환

그림 2.3.2　행렬 \(A\)에 의한 \(q_0\)의 선형 변환

변환된 벡터를 원래 벡터와 함께 생성하는 공간 안에서 직교화한다. 즉, 그림 2.3.3과 같이 Gram-Schmidt 직교화를 수행한다. 그 결과 얻은 벡터를 \(r_1\)이라 하고, 이를 정규화하여 \(q_1\)을 얻는다.

\(q_0\)에 직교하는 벡터 \(q_1\)

그림 2.3.3　\(q_0\)에 직교하는 벡터 \(q_1\)

같은 알고리즘으로 \(q_1\)에서 \(q_2\)를 구한다. 이때 \(q_2\)는 \(q_1\)과 \(q_0\) 모두에 직교한다(그림 2.3.4).

\(q_1\)과 \(q_0\)에 직교하는 벡터 \(q_2\)

그림 2.3.4　\(q_1\)과 \(q_0\)에 직교하는 벡터 \(q_2\)

이와 같이 Lanczos 법은 Krylov 부분공간의 정규직교 기저를 차례로 구성한다. 이론적으로는 대상 고유값 문제의 대칭성을 이용하여 이 반복을 최근의 기저 벡터를 사용하는 삼항 점화식으로 나타낼 수 있다.

한편 FrontISTR 구현에서는 유한 정밀도 연산으로 인한 직교성 손실을 억제하기 위해 이미 구한 Lanczos 기저 벡터에 대해 \(M\)-내적을 이용한 재직교화를 수행한다.

삼중대각화

FrontISTR의 Lanczos 반복에서는 앞 절에서 설명한 \(M\)-내적에 대해 기저 벡터를 정규직교화하므로

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

위 관계가 성립한다.

식 \(\eqref{eq:2.3.10}\)의 행렬 \(A\)를 사용하면 이론적인 Lanczos 반복은

\[\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}\)을 이미 구한 Lanczos 기저 벡터에 대해 \(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}\]

위 관계가 성립한다.

Lanczos 반복으로 얻은 \(m\)개의 기저 벡터를 다음과 같이 모으면

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

유한 횟수의 Lanczos 반복 후에는

\[\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\)성분만 1인 \(m\)차원 단위 벡터이며

\[\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\)의 고유값을 계산하면 원래의 대규모 고유값 문제의 고유값을 근사할 수 있다.

관련 항목