固有値解析

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

質量モデル

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

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

両方式の公称並進総動的質量は \( \sum W_{\text{mass},i}/g \) で一致します。動的質量は物理質量（鋼材は物理質量密度 7.85 t/m³）に基づき、地震用重量 \( W \)（設計重量。鋼材は 78.5 kN/m³）とは一致しません。

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

質点系（串団子）の層質量は、この質量行列とは別に、階生成が物理質量相当重量の全量から Story::dynamic_mass として保存した値を用います（質点系解析の「質量」）。

記号・単位

算定式

要求モード数 \( p_{\text{req}} \) は縮約後自由度数 \( n \) で切り詰め、解析対象のモード数 \( p \) とします:

\[ p = \min(p_{\text{req}}, n) \]

部分空間サイズ \( q \) は次で算定します:

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

上式は \( p_{\text{req}} \geq 1 \) に適用します。要求モード数 \( p_{\text{req}} = 0 \) のときは上式によらず \( q = 0 \) とし、部分空間を構成しません（固有値解析は行わず、空の結果を返します）。

この式は Bathe 2013 式(19) の \( q = \max(p+8, 2p) \) に、縮約後自由度数 \( n \) の上限を課したものです。部分空間反復は本来、Sturm 列チェックによる固有振動数の欠落検出と部分空間の自動拡大を併用しますが、本実装では Sturm 列チェックを行わないため、小さい部分空間でモードを取りこぼすのを避けて \( 2p \) と \( p+8 \) の大きい方（\( \max \) 側）を採ります。

要求モード数 \( p_{\text{req}} \) が縮約後自由度数 \( n \) を超える場合は、自由度数まで切り詰めて \( p = n \) として実行します（返るモード数は要求モード数を超えません）。一方、質量を持つ方向の数（質量ランク）が解析対象のモード数 \( p \) に満たない場合は、要求モード数を切り詰めるのではなくエラーとして通知します（質量が有効な独立自由度が不足しているため）。

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

\[ 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}}) \)、合同変換のため固有値は不変）することで、並進質量と回転慣性のように単位スケールが大きく異なる質量が混在しても、質量の有無の判定が単位系に依存しません。質量を持たない方向（対角化後の質量対角成分が相対許容誤差未満の方向）には固有振動数が存在しないため解の対象から除外し、質量を持つ方向の数（質量ランク）が解析対象のモード数 \( p \) に満たない場合はエラーとして通知します。

収束は \( |\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 分解（線形代数ソルバ）を（部分空間サイズ×反復回数）回の前進・後退代入で再利用します。固有値解析のソルバは、自由度規模によらず常に疎 Cholesky 直接法です。部分空間反復は同一の分解を多数回の求解で再利用する構造のため、求解のたびに反復計算をやり直す反復法（PCG）は用いません（線形代数ソルバの AUTO 選択の例外）。静解析と続けて実行する場合、静解析の準備が直接法で分解している規模（自由度 50,000 未満）ではその分解をそのまま使い回し、固有値解析のための再分解は行いません。静解析側が反復法（PCG）を選ぶ規模（自由度 50,000 以上）では、固有値解析専用に直接法で 1 回だけ分解します。

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

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

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