Files
김경종 bd50e09e36
Tests / Hermetic test suite (push) Has been cancelled
Tests / Skill frontmatter validation (push) Has been cancelled
add documents
2026-06-02 11:38:52 +09:00

271 lines
15 KiB
Markdown
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<!-- source-page: 361 -->
# 2-2 Lanczos 반복법
Lanczos반복법은 Krylov 부분공간(subspace) span( $V_1, V_2, ..., V_k$ ) 을 생성하는 과정을 통하여 발생하는 삼중 대각행렬 $T_k$ 를 이용하여 고유치의 근사값을 구하는 방법이다. 구조물의 진동해석에서 발생하는 고유치 문제 식 (2.1.7)에 Lanczos반복법을 효과적으로 적용하기 위해서는 고유치 $\lambda_m$ 을 $\lambda_m = \sigma + 1/\theta_m$ 로 치환하여야 한다. 이러한 방법을 Shift-invert 기법이라 하며, $\sigma$ 는 첫번째 고유치로 예상되는 값이다. Shift-invert 기법을 적용한 Lanczos 반복계산 과정을 간단하게 요약하면 다음과 같다.
- 첫 번째 반복계산의 경우 블록(block) 벡터의 초기값 $V_{1}$ 을 가정
- 질량 행렬 곱셈 $U_{k} = MV_{k}$
- 선형 연립 방정식 풀이 $(\mathbf{K}-\sigma\mathbf{M})\mathbf{W}_{k}=\mathbf{U}_{k}$
- 직교화 $\mathbf{W}_{k}^{*} = \mathbf{W}_{k} - \mathbf{V}_{k-1} \mathbf{B}_{k-1}^{T}$
- 행렬 $C_{k}$ 계산 $C_{k}=V_{k}MW_{k}^{*}$
- 직교화 $\mathbf{W}_{k}^{**} = \mathbf{W}_{k}^{*} - \mathbf{V}_{k} \mathbf{C}_{k}$
- 블록 벡터 내 정규화 $\mathbf{W}_{k}^{**} = \mathbf{V}_{k+1} \mathbf{B}_{k}$
midas FEA에서는 효과적인 고유치 계산을 위하여 블록 벡터 $V_{m}$ 를 이용하며, 이를 블록 Lanczos 방법이라 한다. 위의 반복계산 과정을 통하여 발생하는 블록 삼중 대각 행렬 $T_{k}$ 는 다음과 같다.
$$
\mathbf {T} _ {k} = \left[ \begin{array}{c c c c c} \mathbf {C} _ {1} & \mathbf {B} _ {1} ^ {T} & & & \\ \mathbf {B} _ {1} & \mathbf {C} _ {2} & \ddots & & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \mathbf {C} _ {k - 1} & \mathbf {B} _ {k - 1} ^ {T} \\ & & & \mathbf {B} _ {k - 1} ^ {T} & \mathbf {C} _ {k} \end{array} \right] \tag {2.2.1}
$$
$T_{k}$ 를 이용하여 고유치 문제 $T_{k}\psi_{m}=\theta_{m}^{*}I$ 를 풀게 되면 $\lambda_{m}^{*}=\sigma+1/\theta_{m}^{*}$ 을 이용하여 $\lambda_{m}^{*}$ 을 구할 수 있으며, $\lambda_{m}^{*}$ 은 본래의 고유치 문제 식 (2.1.7)의 근사값이다. $V_{k}$ 의 블록 크기를 $N_{b}$ 이라 할 때, 반복계산 과정의 반복 회수가 늘어날 수록 $T_{k}$ 의 크
<!-- source-page: 362 -->
기는 Nb 만큼 커지게 되며, \* λm 는 λm 으로 수렴한다. 고유모드 φ m 의 근사값 \* φm 는다음 식을 통하여 구할 수 있으며, \* λ 과 함께 수렴한다.
$$
\phi_ {m} ^ {*} = \left[ \begin{array}{l l l l l} \mathbf {V} _ {1} & \mathbf {V} _ {2} & \mathbf {V} _ {3} & \dots & \mathbf {V} _ {k} \end{array} \right] \psi_ {m} ^ {*} \tag {2.2.2}
$$
Lanczos 반복법으로 계산된 고유치와 고유모드의 수렴 여부는 다음과 같이 확인할 수 있다.
$$
\frac {\left\| \mathbf {K} \phi_ {m} ^ {*} - \lambda_ {m} ^ {*} \mathbf {M} \phi_ {m} ^ {*} \right\|}{\left\| \mathbf {K} \right\|} \leq \varepsilon \tag {2.2.3}
$$
여기서, . 은 2-norm을 의미한다. midas FEA에서는 16 ε 2.22 10 = × 를 이용한다.
<!-- source-page: 363 -->
# 2-3 Subspace 반복법
Subspace 반복법은 부분공간 $E_{k}$ 를 이루는 $N_{s}$ 개의 벡터 $X_{k}$ 가 반복계산과정을 통하여 고유모드 $\left[\phi_{1} \quad \phi_{2} \quad \ldots \quad \phi_{N_{s}}\right]$ 에 수렴하도록 하는 방법이다. 부분공간 $E_{k}$ 를 구성하는 벡터 $X_{k}$ 를 반복적으로 계산해 나가는 과정은 다음과 같다.
\- 첫 번째 반복계산의 경우 $E_{1}$ 을 구성하는 $N_{s}$ 개의 초기 벡터 $X_{1}$
-선형 연립 방정식 풀이 $\mathbf{K}\mathbf{Y}_{k}=\mathbf{M}\mathbf{X}_{k}$
\- 강성 행렬의 투영(projection) $\mathbf{K}_{k+1}=\mathbf{Y}_{k}^{T}\mathbf{K}\mathbf{Y}_{k}$
\- 질량 행렬의 투영(projection) $\mathbf{M}_{k+1}=\mathbf{Y}_{k}^{T}\mathbf{M}\mathbf{Y}_{k}$
\- 투영된 고유치 문제 풀이 $K_{k+1}Q_{k+1}=M_{k+1}Q_{k+1}\Lambda_{k+1}$
$$
\mathbf {Q} _ {k + 1} = \left[ \begin{array}{c c c c} \psi_ {1} ^ {*} & \psi_ {2} ^ {*} & \dots & \psi_ {N _ {s}} ^ {*} \end{array} \right]
$$
$$
\boldsymbol {\Lambda} _ {k + 1} = \left[ \begin{array}{c c c c} \lambda_ {1} ^ {*} & & & \\ & \lambda_ {2} ^ {*} & & \\ & & \ddots & \\ & & & \lambda_ {N _ {s}} ^ {*} \end{array} \right]
$$
$$
- \mathbf {X} _ {k + 1} \text { 의 계산 } \quad \mathbf {X} _ {k + 1} = \mathbf {Y} _ {k} \mathbf {Q} _ {k + 1}
$$
위의 반복계산과정을 통하여 계산되는 $\lambda_{n}^{*}$ 와 $X_{k}$ 는 다음과 같이 고유치와 고유모드로 각각 수렴하게 된다.
$$
\lambda_ {m} ^ {*} \rightarrow \lambda_ {m}, \quad \mathbf {X} _ {k} \rightarrow \left[\begin{array}{c c c c}\phi_ {1}&\phi_ {2}&\dots&\phi_ {N _ {s}}\end{array}\right] \tag {2.3.1}
$$
$X_{k}$ 가 $N_{s}$ 개의 벡터로 이루어져 있을 때, 반복계산과정의 반복 회수와 관계 없이 $K_{k}$ 와 $M_{k}$ 의 크기는 $N_{s} \times N_{s}$ 로 고정된다. Subspace 반복법으로 계산된 고유치의 수렴 여부는 근사값 $\lambda_{m}^{*}$ 이 변화하는 정도에 의하여 다음과 같이 판단한다.
<!-- source-page: 364 -->
Part 4 Linear Analysis
$$
\left\| \frac {\lambda_ {m} ^ {(k + 1)} - \lambda_ {m} ^ {(k)}}{\lambda_ {m} ^ {(k + 1)}} \right\| \leq \varepsilon \tag {2.3.2}
$$
<!-- source-page: 365 -->
# 2-4 관련기능
# 2-4-1 질량 행렬
고유치 해석에서 사용할 수 있는 질량행렬의 종류에는 분포(consistent) 질량과 집중(lumped) 질량이 있다. 분포 질량을 사용하게 되면 이론 해에 비하여 고유치가크게 계산되며, 집중 질량의 경우는 일반적으로 고유치가 작게 계산된다. 분포질량은 집중질량에 비하여 요소수에 따른 고유치의 수렴성이 더 좋은 것으로 알려져있으나, 계산량과 메모리 요구량이 많은 단점이 있다. 이와 같은 특징들은 그림2.4.1을 통하여 알 수 있다.
![](images/page-365_2151cb534eec83186770c6f7025ac8f727cf11cb295d4d8f0888221e23fa0354.jpg)
<details>
<summary>line</summary>
| Number of Element | Lumped Mass | Consistent Mass |
| ----------------- | ----------- | --------------- |
| 0 | Low | High |
| 1 | Medium | Medium |
| 2 | High | Low |
| 3 | High | Low |
| 4 | High | Low |
| 5 | High | Low |
| 6 | High | Low |
| 7 | High | Low |
| 8 | High | Low |
| 9 | High | Low |
| 10 | High | Low |
</details>
그림 2.4.1 집중질량과 분포질량의 진동수 비교
<!-- source-page: 366 -->
# 2-4-2 Lanczos 반복법 관련기능
midas FEA에서 사용하는 Lanczos 반복법은 두 개의 진동수 1f , 2f 에 의해 계산하고자 하는 진동수의 관심영역(range of interest)을 결정할 수 있다. 그리고 진동수의 관심 영역은 Shift-invert 기법에서 고유치 예상값을 $\sigma = ( 2 \pi f _ { 1 } ) ^ { 2 }$ 로 설정함으로써 계산에 반영할 수 있다. 그림 2.4.2는 1f 과 2f 의 관계에 따라 계산하게 되는고유치의 영역과 순서를 그린 것이다. 계산하고자 하는 고유치 개수를 $N _ { f }$ 이라 할때, 1f , 2f 의 크기에 따라 다음과 같이 고유치를 계산한다.
- 1 2f = f 인 경우 : 1f 에 가까운 진동수 $N _ { f }$ 개를 계산
- 1 2f < f 인 경우 : 영역 1 2[ , ] f f 에서 1f 에 가까운 진동수 $N _ { f }$ 개를 계산
- 1 2f > f 인 경우 : 영역 2 1[ , ] f f 에서 1f 에 가까운 진동수 $N _ { f }$ 개를 계산
![](images/page-366_f7cc228ecf2ac917ddf6037ba874c30bb1dc477738c0bdcd512cb6faba515316.jpg)
<details>
<summary>text_image</summary>
f₁=f₂
④ ①②③
←→
f
</details>
![](images/page-366_b00ef754b887d00ff6247caade26a7e2c9f4ee27399bf6aa4b55eb69353c7449.jpg)
<details>
<summary>text_image</summary>
f₁
①②
③④
f₂
f
</details>
![](images/page-366_d2c28fe27574e31158106430a02f4a9d10b1bdcba1e694f166a4de4e8f5448bd.jpg)
<details>
<summary>text_image</summary>
f₂
③②①
f₁
f
</details>
그림 2.4.2 고유치 검색 방향
진동수의 관심영역 1f , 2f 의 기본값은 f f 1 2= = 0 이다.
<!-- source-page: 367 -->
# 2-4-2-1 Sturm Sequence
반복법을 이용하여 고유치를 계산하게 되면 낮은 차수의 고유치보다 높은 차수의고유치가 먼저 수렴하는 경우가 발생한다. midas FEA에서는 이와 같은 현상을 고려하여 진동수의 누락을 방지하는 기능을 제공하며, 그 방법은 다음과 같다.Lanczos 반복계산 과정 중에는 다음과 같은 연립 방정식의 해법이 필요하다.
$$
(\mathbf {K} - \sigma \mathbf {M}) \mathbf {W} _ {k} = \mathbf {U} _ {k} \tag {2.4.1}
$$
이 방정식의 풀이를 위하여 T LDL 행렬 분해를 수행하고, 행렬 D의 대각항 중 음수의 개수를 계산한다. 이 개수가 σ 보다 작은 고유치 λn 의 개수이다. 예를 들어$\boldsymbol { \sigma } = \lambda _ { N _ { f } } + \delta \underline { { \boldsymbol { \circ } } } \underline { { \boldsymbol { \mathbf { \Pi } } } }$ 두고 행렬을 분해하면 $\lambda _ { N _ { f } } + \delta$ 보다 작은 고유치 개수를 알 수 있다. $\lambda _ { N _ { f } } + \delta$ 가 보다 작은 고유치 개수가 $N _ { f }$ 보다 크면, 누락된 고유치가 모두 수렴할 때까지 Lanczos 반복계산 과정을 지속한다. 이 때 δ는 다음의 관계를 만족하도록 결정한다.
$$
\lambda_ {N _ {f}} < \lambda_ {N _ {f}} + \delta < \lambda_ {N _ {f} + 1} \tag {2.4.2}
$$
midas FEA에서는 “Sturm Sequence Check”를 수행할 경우 Lanczos 블록 벡터$\mathbf { V } _ { k }$ 의 크기를 $N _ { b } = 7$ 로 하여 계산 속도를 높이도록 하였다. 그러나 T LDL 분해를추가적으로 수행해야 하는 계산 부담을 고려하여 “Sturm Sequence Check”를 생략할 수 있도록 선택사항을 제공하고 있으며, 이 경우에는 $N _ { b } = 1$ 로 계산하기 때문에 고유치의 누락을 최소화할 수 있다.
# 2-4-2-2 강체 모드
진동 특성을 분석하고자 하는 구조물의 경계조건이 부족한 경우에는 K 가 특이(singular) 행렬이 되기 때문에 σ = 0 으로 계산을 진행할 수 없다. 이러한 경우에는 0보다 작은 σ 를 가정하여 Shift-invert 기법을 적용하며, 계산된 고유치는 부족한 경계조건 개수 만큼의 0을 포함한다. midas FEA에서는 행렬 K 와 M 의 대각항을 이용하여 1차 모드로 예상되는 진동수 $\overline { { \lambda } } _ { 1 } \equivq$ 예측하고, $\sigma = - \overline { { \lambda } } _ { 1 }$ 로 설정한다.
<!-- source-page: 368 -->
# 2-4-3 Subspace 반복법 관련기능
midas FEA에서 사용하는 Subspace 반복법에서는 Subspace $\mathbf { X } _ { k } \ \underline { { \circ } } | \ \exists 7 | \ N _ { s }$ ,최대 반복 회수 $N _ { I }$ , 그리고 식 (2.3.2)의 수렴 허용오차 $\varepsilon \ { \stackrel { \circ } { \equiv } }$ 설정할 수 있다.$N _ { s }$ 의 기본값은 0이며, 실제로 계산에 사용되는 Subspace의 크기 $ { N _ { s } } \equiv \ { \mathbb { L } } | \mathbf { { \frac { \mathrm { ~ o ~ } } { \mathrm { ~ \scriptstyle { \Xi ~ } } } } } { \mathbf { \vec { \imath } } } |$ 같다.
$$
N _ {s} = \max \left\{N _ {s}, \min (2 N _ {f}, N _ {f} + 8) \right\} \tag {2.4.3}
$$
Subspace 반복계산 과정의 최대 반복 회수 기본값은 $N _ { \mathit { I } } = 3 0 \ 0 | \mathbb { H }$ , 수렴 허용오차 기본값은 $\varepsilon = 1 . 0 \times 1 0 ^ { - 6 } 0$ 다.
<!-- source-page: 369 -->
# Chapter 3. Time History Analysis
# 3-1 개요
시간이력해석은 구조물에 동적하중이 작용할 때 식(3.1.1)과 같은 동적 평형방정식의 해를 구하는 것이다.
$$
\mathbf {M} \ddot {\mathbf {u}} (t) + \mathbf {C} \dot {\mathbf {u}} (t) + \mathbf {K} \mathbf {u} (t) = \mathbf {F} (t) \tag {3.1.1}
$$
여기서,
M : 질량행렬
C : 감쇠행렬
K : 강성행렬
$\mathbf{F}(t)$ : 동적하중
u(t) : 상대변위
$\dot{\mathbf{u}}(t)$ : 속도
$\ddot{\mathbf{u}}(t)$ : 가속도
식 (3.1.1)에서 $\mathbf{F}(t)=0$ 이면 자유진동에 대한 방정식이 되고, C=0 인 조건을 추가하면 비감쇠 자유진동방정식이 된다. 그리고 $\mathbf{F}(t)$ 가 임의 시간에 대한 기진력(또는 기진변위, 속도, 가속도 등)이면 강제진동 해석문제가 된다. midas FEA에서는 식(3.1.1)로 표현되는 시간이력해석을 위해 모드중첩법(modal superposition method)과 직접적분법(direct integration method)을 사용한다
<!-- source-page: 370 -->
# 3-2 모드중첩법
모드종첩법에서는 고유모드의 직교성을 이용하여, 식(3.1.1)을 모드별로 서로 독립된 방정식으로 분해하여 계산을 수행한다. 이를 위하여 모드종첩법에서는 감쇠행렬은 질량행렬과 강성행렬의 선형조합으로 이루어진다는 가정을 사용한다. 그리고 모드별 응답이 얻어지면 이들을 식(3.2.1)와 같이 선형조합하여 해를 구한다.
$$
\mathbf {u} (t) = \boldsymbol {\Phi} \mathbf {Y} (t) = \sum_ {i = j} ^ {n} \boldsymbol {\phi} _ {i} Y _ {i} (t) \tag {3.2.1}
$$
식(3.2.1)을 식(3.1.1)에 대입하면 식(3.2.2)와 같이 나타낼 수 있다.
$$
\mathbf {M} \boldsymbol {\Phi} \ddot {\mathbf {Y}} (t) + \mathbf {C} \boldsymbol {\Phi} \dot {\mathbf {Y}} (t) + \mathbf {K} \boldsymbol {\Phi} \mathbf {Y} (t) = \mathbf {F} (t) \tag {3.2.2}
$$
식(3.2.2)에 $\phi_{m}^{T}$ (m번째 모드형상)을 곱하면 식(3.2.3)가 얻어진다.
$$
\boldsymbol {\phi} _ {m} ^ {T} \mathbf {M} \boldsymbol {\Phi} \ddot {\mathbf {Y}} (t) + \boldsymbol {\phi} _ {m} ^ {T} \mathbf {C} \boldsymbol {\Phi} \dot {\mathbf {Y}} (t) + \boldsymbol {\phi} _ {m} ^ {T} \mathbf {K} \boldsymbol {\Phi} \mathbf {Y} (t) = \boldsymbol {\phi} _ {m} ^ {T} \mathbf {F} (t) \tag {3.2.3}
$$
고유모드의 질량행렬과 강성행렬에 대한 직교성은 다음과 같이 표현할 수 있다.
$$
\phi_ {i} ^ {T} \mathbf {M} \phi_ {j} = m _ {i j} \delta_ {i j} \tag {3.2.4}
$$
$$
\boldsymbol {\phi} _ {i} ^ {T} \mathbf {K} \boldsymbol {\phi} _ {j} = k _ {i j} \boldsymbol {\delta} _ {i j}
$$
여기서,
$$
\delta_ {i j} \quad : \text { 크로네커 델타(Kronecker delta) } \quad (i = j; \delta_ {i j} = 1, i \neq j; \delta_ {i j} = 0)
$$
고유모드의 질량 및 강성행렬에 대한 직교성을 이용하면, 동적평형방정식(3.2.2)은식(3.2.5)과 같이 서로 독립된 n개의 독립된 방정식으로 분해할 수 있다.
$$
\boldsymbol {\phi} _ {m} ^ {T} \mathbf {M} \boldsymbol {\phi} _ {m} \ddot {Y} (t) + \boldsymbol {\phi} _ {m} ^ {T} \mathbf {C} \boldsymbol {\phi} _ {m} \dot {Y} (t) + \boldsymbol {\phi} _ {m} ^ {T} \mathbf {K} \boldsymbol {\phi} _ {m} Y (t) = \boldsymbol {\phi} _ {m} ^ {T} \mathbf {F} (t) \quad (m = 1, 2, \dots \dots , n) \tag {3.2.5}
$$
따라서 모드별로 분해된 m개의 식은 일반화 좌표계(general coordinate)에서의 단자유도계에 대한 동적평형방정식이 된다. 그리고 식(3.2.5)를 변형하면 식(3.2.6)