5.3 固有値解析

固有周期とモードを、部分空間反復法（Bathe）と一般化 Jacobi 法による射影固有値問題の求解で求めます。

質量モデル

全体質量行列 \( M \) の質量源は「部材密度による整合質量（Consistent）」と「節点集中質量（剛床マスターの質点等）」の 2 つで、その組合せはモデルの質量方式に従います:

• 補正質点方式（既定）: 部材密度による分布質量に加え、剛床マスター節点へ「地震用重量のうち分布質量として計上されない分」（床・仕上げ・積載・二次部材・雑壁など）を補正質点として与えます。階の生成がマスター質量を \( m = (W - W_{\text{frame}})/g \)（水平 2 方向）、回転慣性を \( J = \sum m_i r_i^2 \)（\( r_i \) は慣性力重心からの平面距離）で算定します。階の合計質量は地震用重量 \( W \) と一致し、部材の分布質量による鉛直・局部の振動モードも保たれます。
• 質点のみ方式: 質量は節点集中質量のみを用い、部材密度による分布質量は算入しません（水平質点系モデル化）。剛床マスターには地震用重量の全量 \( W/g \) と回転慣性を与えるため、鉛直方向・局部の振動モードは表現されません。

質量方式は固有値解析・時刻歴応答解析・精算固有周期（Ai 分布の T 算定）で共通に適用されます。

算定式

部分空間サイズ（\( p = \text{モード数} \)）:

\[ q = \min(2p, p+8) \]

反復は次の予測から射影を作り、射影一般化固有値問題を解いて更新します:

\[ K y = M x \]

\[ \bar{K} = y^\mathsf{T} K y, \bar{M} = y^\mathsf{T} M y \]

\[ x_{\text{new}} = y \cdot V \]

射影一般化固有値問題 \( \bar{K} V = \bar{M} V \Theta \) は、一般化 Jacobi 法（Bathe）により \( \bar{K} \)・\( \bar{M} \) を同時対角化して解きます。半正定値の \( \bar{M} \)（部分空間サイズが質量ランクを超える場合に必ず生じる）でも正則な合同変換のみで進むため、質量を持たない方向が混在しても数値的に安定します。求解に先立ち、両行列を \( \bar{K} \) の対角で対称スケーリング（\( S = \operatorname{diag}(1/\sqrt{\bar{K}_{ii}}) \)、合同変換のため固有値は不変）することで、並進質量と回転慣性のように単位スケールが大きく異なる質量が混在しても、質量の有無の判定が単位系に依存しません。質量を持たない方向（対角化後の質量対角成分が相対許容誤差未満の方向）には固有振動数が存在しないため解の対象から除外し、質量を持つ方向の数（質量ランク）が要求モード数に満たない場合はエラーとして通知します。

収束は \( |\theta - \theta_{\text{prev}}| < 10^{-10} \cdot \max(\theta, 1) \) で判定します。周期および固有ベクトルの M 正規化は次のとおりとします:

\[ \text{period} = 2\pi/\sqrt{\omega^2} \]

\[ \varphi^\mathsf{T} M \varphi = 1 \]

固有ベクトル \( \varphi \) は拘束縮約後の独立自由度空間で求まるため、モード形状の表示・出力にあたっては縮約を逆変換して全節点の 6 成分へ展開します。剛床などの従属自由度には、マスター自由度の剛体変位関係（面内並進＋回転）に従った値が入ります。

参加係数および有効質量は次で算定します:

\[ \beta = \frac{\varphi^\mathsf{T} M r}{\varphi^\mathsf{T} M \varphi} \]

\[ M_{\text{eff}} = \frac{(\varphi^\mathsf{T} M r)^2}{\varphi^\mathsf{T} M \varphi} \]

疎行列演算と計算コスト

全体剛性行列 \( K \)・全体質量行列 \( M \) は、立体骨組ではほとんどの成分がゼロの疎行列になります（各自由度は同じ節点・隣接節点の自由度としか連成しないため）。固有値解析の反復に現れる行列演算は、この疎行列の非ゼロ要素だけを用いて計算します:

1. 射影行列の算定：部分空間への射影 \( \bar{K} = Y^\mathsf{T} K Y \)、\( \bar{M} = Y^\mathsf{T} M Y \) は、まず疎行列ベクトル積 \( Z = K Y \)（列ごとに非ゼロ要素のみを走査。コストは非ゼロ数 \( \mathrm{nnz} \) に比例）を求め、続いて \( Y^\mathsf{T} Z \) を計算する 2 段構成とします。1 反復あたりのコストは \( O(\mathrm{nnz} \cdot q + n q^2) \)（\( n \): 自由度数、\( q \): 部分空間サイズ）で、モデル規模に対しほぼ線形にしか増えません。
2. 固有ベクトルの M 正規化・参加係数：\( \varphi^\mathsf{T} M \varphi \) や \( \varphi^\mathsf{T} M r \) も同様に疎行列ベクトル積で算定します。
3. 分解の再利用：反復中の \( y = K^{-1} M x \) は、解析開始時に 1 回だけ行う剛性行列の疎 Cholesky 分解（5.1 線形代数ソルバ）を（部分空間サイズ×反復回数）回の前進・後退代入で再利用します。固有値解析のソルバは、自由度規模によらず常に疎 Cholesky 直接法です。部分空間反復は同一の分解を多数回の求解で再利用する構造のため、求解のたびに反復計算をやり直す反復法（PCG）は用いません（5.1 の AUTO 選択の例外）。静解析と続けて実行する場合、静解析の準備が直接法で分解している規模（自由度 50,000 未満）ではその分解をそのまま使い回し、固有値解析のための再分解は行いません。静解析側が反復法（PCG）を選ぶ規模（自由度 50,000 以上）では、固有値解析専用に直接法で 1 回だけ分解します。

このため固有値解析の所要時間は、モデル規模（自由度数と非ゼロ数）にほぼ比例して増える程度に収まり、小規模モデルでは即座に完了します。

地震荷重（Ai 分布）との連携

地震荷重の設計用一次固有周期 T の算定法（1.1 地震力）で「精算（固有値解析）」を選択した場合は、本解析で求めた 1 次モードの固有周期を用います。固有周期は加力方向（X/Y）によらず同一のため、算定は 1 回で足ります。精算周期を使うには、解析タブで固有値解析を先に実行しておく必要があります（未実行の場合は画面に案内が表示されます）。既定の「略算」は告示の略算式 T = h(0.02+0.01α) を用いるため、固有値解析の実行は不要です。

実装：dynamic::eigen（eigen/mod.rs::solve_eigen / solve_eigen_with_solver）が算定します。 EIGEN_TOL = 1e-10、EIGEN_MAX_ITER = 200 とし、部材の分布質量は整合質量（Consistent）で組み立てます（算入の要否は上記の質量方式に従います）。整合質量行列は要素局所系で構成し、剛性行列と同様に全体系へ変換（\( M_{\text{global}} = R^\mathsf{T} M_{\text{local}} R \)）してから全体自由度へ組み込みます（軸方向と曲げ方向で係数が異なり回転不変ではないため。部材の向きに依らず同一の固有周期が得られます）。