5.5 時刻歴応答解析

支配方程式（基盤一様加振、相対変位形式）を解きます。全体質量行列 \( M \) は固有値解析と同一の質量モデル（5.3 固有値解析の質量モデル: 部材の整合質量と剛床マスター等の節点集中質量を、モデルの質量方式に従って組み合わせる）で組み立てます。

\[ M \cdot \ddot{u} + C \cdot \dot{u} + K \cdot u = -M \cdot r \cdot \ddot{x}_g(t) \]

5.5.1 線形時刻歴（Newmark-β）

線形の時刻歴応答を Newmark-β 法で解きます。 パラメータは平均加速度法（\( \beta = 1/4, \gamma = 1/2 \)、無条件安定）と線形加速度法（\( \beta = 1/6, \gamma = 1/2 \)、条件付安定）を取り、解析は常に平均加速度法で実行します。 時間刻みによらず安定なため、波形の刻みをそのまま使えるためです。 時間刻み \( \Delta t \) は、解析設定で指定した値が正であればそれを、指定がなければ入力波形のサンプリング間隔を使います。

算定式

有効剛性（\( c_1 = 1/(\beta \Delta t^2), c_2 = \gamma/(\beta \Delta t) \)）:

\[ \hat{K} = K + c_2 \cdot C + c_1 \cdot M \]

次を解き、加速度と速度を更新します:

\[ u_{n+1} = \hat{K}^{-1} \cdot p_{\text{eff}} \]

\[ a_{n+1} = c_1(u_{n+1} - u) - c_3 \cdot v - c_4 \cdot a \]

\[ v_{n+1} = v + \Delta t((1-\gamma)a + \gamma \cdot a_{n+1}) \]

初期加速度は次から求めます:

\[ M \cdot a_0 = -C \cdot v_0 - K \cdot u_0 - p(0) \]

実装：dynamic::timehistory（timehistory/mod.rs::linear_time_history_with_state）が算定します。

5.5.2 非線形時刻歴（Newmark-β ＋ Newton-Raphson）

非線形の時刻歴応答を、Newmark-β と Newton-Raphson 反復（commit/rollback）で解きます（構造動力学の標準的な非線形時刻歴解法）。 積分パラメータと時間刻みの決め方は 5.5.1 と同じで、平均加速度法（\( \beta = 1/4, \gamma = 1/2 \)）を使います。

算定式

各反復で接線剛性 Kt と接線比例減衰を再構成し、有効剛性を次とします:

\[ K_{\text{eff}} = K_t + c_2 \cdot C + c_1 \cdot M \]

残差を次とします（\( f_{\text{int}} \) は要素の復元力に加え、支点ばね（Node::support_spring） の内力 \( k_i u_i \) も含みます）:

\[ r = p - f_{\text{int}} - C \cdot v - M \cdot a \]

収束は \( \|r\| < \text{tol} \cdot r_{\text{ref}} \) で判定します。基準ノルム \( r_{\text{ref}} \) は 動的釣り合いの各項のノルムの最大とし、解析中の最大値に対する下限を設けます。

\( p_{\text{dyn}} \) は外力 \( p \) から長期荷重成分 \( f_0 \) を除いた地震動由来の 動的外力（\( p = p_{\text{dyn}} + f_0 \)）です。\( f_0 \) を含めないのは、長期荷重が 卓越するモデルで収束判定が過度に緩まないようにするためです。

慣性力・減衰力も基準に含めるのは、地動加速度がゼロを横切る時刻で基準が消えるのを 防ぐためです。動的外力だけを基準にすると、その時刻で分母が最低値の 1 まで落ち、判定が 「残差が tol（N）を下回ること」という絶対値判定に化けます。実建物では内力・慣性力が \( 10^7 \) N 規模あります。有効剛性の線形解が持つ残差は、その条件数と計算機イプシロンに 比例して \( 10^{-1} \) N 規模になります。つまり到達できない閾値です。正弦波なら 毎周期 2 回通る時刻であり、実地震波でも頻繁に起きます。地動が 0 でも構造は動いており 慣性力・減衰力は大きいので、3 項の最大を採れば基準は消えません。

3 項の最大だけでは、応答が減衰しきった時刻で同じ破綻が戻ります。長い波形の末尾では 3 項すべてが 1 N を下回り、また床まで落ちるためです。そこで解析中に観測した力のスケールの 最大値に対する下限を設け、その 1% を切った領域では基準をそれ以上下げません。

\[ r_{\text{ref}} = \max(\|p_{\text{dyn}}\|,\ \|M \cdot a\|,\ \|C \cdot v\|,\ 0.01 \cdot r_{\text{ref,max}},\ 1) \]

応答が最大の 1% を切った領域は工学的に無視できる大きさで、そこで判定を続けても意味が ありません。実測では、最大 6.1e7 N の解析の末尾で力が 6e-4 N まで落ちたところ、履歴則が 微小振幅でチャタリングして残差が 1e-4 N から下がらなくなりました。1% はこれを収束と 認めるのに必要な余裕から定めています。応答が大きい時刻では瞬間値のほうが大きいため、 この下限は効きません。

最低値の 1 は、解析を通して一度も動いていないときに 0 除算を避けるための床です。 反復の最大回数と許容誤差 tol は解析タブ「時刻歴応答」の「Newton反復」欄で指定します （既定は最大 50 回・tol \( 10^{-6} \)）。

次を解いて要素状態を更新します:

\[ \delta u = K_{\text{eff}}^{-1} \cdot r \]

反復が上限に達しても収束しなかったステップは、その時点の試行状態で確定して解析を続行し、 非収束ステップ数を結果画面に注記します。途中まで解けた応答を捨てるより、参考値として 示したうえで信頼性の低下を伝えるほうが判断材料になるためです。この場合の応答値は残差の 収束を確認できていない参考値のため、時間刻み dt を小さくするなどの見直しを検討してください。 質点系（5.5.3）も同じ扱いです。

ただし応答が発散して変位が有限値でなくなった場合は、ステップ開始状態へ戻して解析を 打ち切ります。以降のステップも結果もすべて無効になり、参考値にもならないためです。

実装：timehistory/mod.rs::nonlinear_time_history_analysis が算定します。 コンクリート履歴は動的解析で原点指向型に切替え、制振要素へ set_time_step(dt) を設定します。 位相差入力（ねじれ加振）にも対応します。 非線形時刻歴の幾何剛性は考慮しません。

入力方向と位相差入力

X・Y に加え、同一波形を両方向へ同時入力する X+Y を選べます。 X+Y は同一波形をそのまま両方向へ入れる簡易仕様で、位相差や方向別の波形指定には対応しません。

位相差入力では、位相遅れ時間 \( t = (L \cdot \sin\theta)/V_s \) を求め、位相遅れ方向の並進波を基準波としてねじれ地動加速度を生成します。 \( L \) は矩形基礎長さ [m]、\( V_s \) はせん断波速度 [m/s]、\( \theta \) は入射角 [°] です（既定は \( L = 20 \)、\( V_s = 200 \)、\( \theta = 30 \)）。

長期荷重の初期化

時刻歴の開始前に、長期系荷重ケース（固定・積載等、LoadCaseKind::is_long_term）の外力を 静的 Newton 反復で載荷し（7. プッシュオーバー解析の長期載荷 フェーズと同じ経路）、その変位・応力状態を時刻歴の初期条件とするオプションです（既定で有効）。 軸力に依存する部材耐力（柱の曲げ降伏 \( M_y \)・せん断降伏 \( Q_y \) の軸力項、ファイバー断面の N–M 相関）を、長期軸力を含んだ状態で評価するために用います。

長期荷重ベクトル \( f_0 \) は時刻歴を通じて一定の外力として毎ステップ加算し、動的な初期変位・ 初期速度（波形読込フォームで指定する initial_disp/initial_vel）は「長期解からの増分」として 扱います。したがって、記録される変位・部材内力は長期変位・長期応力を含む全量です。長期系荷重 ケースがないモデルでは何も行いません。このオプションは非線形時刻歴のみに適用され、線形時刻歴 （5.5.1）には適用されません（地震動による応答成分のみを扱います）。

5.5.3 質点系（串団子）非線形時刻歴

質点系（串団子）モデルの非線形時刻歴を解きます。 モデルの組み立て（2 次元／3 次元、線形の層剛性、非線形骨格）は 5.10 質点系解析 です。 2 次元の非線形時刻歴は、プッシュオーバー Q–δ を等包絡面積則でトリリニアに縮約し、せん断型串団子で Newmark 平均加速度と Newton を用い、三重対角 Thomas 法で解きます（質点系縮約の標準的な実務手法）。 3 次元では各階 \( (U_x, U_y, \theta_z) \) の密な接線剛性を Newton で解きます。ねじりばねは線形です。 層間・層せん断・塑性率は X・Y・45°・水平合成の 4 欄です（定義は 5.10）。 減衰・サンプル波の dt／継続／周期／振幅は立体時刻歴とは独立です（既定は 5.10）。

算定式

復元力は最大点指向型（Clough 系）トリリニアとし、等包絡面積則で次を解いて折点を決定します:

\[ A_{\text{tri}}(\delta_2) = A_{\text{actual}} \]

骨格の元になる層 Q–δ 曲線は、プッシュオーバーの性能曲線から抽出します。長期荷重のみを 初期載荷した状態（水平力ゼロ。7. プッシュオーバー解析の 長期載荷フェーズ）の (δ, Q) を曲線の原点として各点から差し引きます。長期荷重による層の 残留水平変形（非対称な荷重配置で生じ得る）を、水平力に対する復元力特性から切り離すためです。

減衰は初期剛性比例 \( a_1 = 2h/\omega_1 \) とします。ω1 は下記「固有値解析」の 1 次モードを 用い、固有値分解が失敗した場合（質量 0 以下の層がある等）のみ、逆反復法による概算へ フォールバックします。

各時刻ステップの Newton 反復は最大 30 回とし、収束を \( \|r\| < 10^{-6} \cdot r_{\text{ref}} \) で判定します （基準ノルム \( r_{\text{ref}} \) の取り方は 5.5.2 と同じです）。反復が上限内に収束しなかった ステップはその時点の試行状態で確定して解析を続行し、結果画面に非収束ステップ数を 警告表示します。この場合の応答値は残差の収束を確認できていない参考値のため、 時間刻み dt を小さくするなどの見直しを検討してください。

実装：squid_n_solver::dynamic::lumped_mass（crates/squid-n-solver/src/dynamic/lumped_mass/）が算定します。

固有値解析

串団子モデル（各層の質量・初期剛性 K1）から、せん断型の一般化固有値問題 \( K x = \omega^2 M x \) を解きます（M は対角、K はせん断型三重対角）。M が対角である ことを利用して標準固有値問題 \( A = M^{-1/2} K M^{-1/2} \) に帰着し、対称行列の固有値 分解（faer の直接解法）で解きます。

• モード数は質点系パネルのモード数設定（既定 3。立体固有値のモード数とは独立。画面上の上限は 30）を用い、 2 次元は層数、3 次元は 3×層数を超える場合はそこまで切り詰めます。
• 2 次元のモード形状は各層の相対変位（下層→上層）で、最上階の値を 1.0 に正規化します。 3 次元は各階 \( (U_x, U_y, \theta_z) \) を頂部水平変位のノルムで正規化します。 立体モデルの固有値解析結果（M 正規化・別の自由度空間）とは正規化基準が異なるため、モード 形状の値どうしは比較せず、周期のみを比較してください。 刺激係数と有効質量は算定しません。
• 3 次元非線形の並進初期剛性は増分骨格の K1 です。静解析の層 Q/δ とは別に、復元力と同じ骨格から組みます。
• 層に質量が 0 以下のものがある場合は算定できません（実体のない極端な高周波モードが 紛れ込むため）。3 次元では回転慣性 \( J \le 0 \) の階も算定できません。

5.5.4 層応答の集計（層せん断力・層せん断力係数・階加速度・階速度・階変位）

5.5.1〜5.5.2 のいずれの経路でも、応答の詳細記録として各時刻の層応答を階（model.stories、下層→上層）ごとに集計します。

層せん断力（慣性力ベース）

節点慣性力ベクトル \( f_{\text{abs}} = M \cdot \ddot{u} + \ddot{x}_g \cdot M \cdot r \) （\( \ddot{u} \) は相対加速度ベクトル、\( \ddot{x}_g \) は地動加速度、\( M \cdot \ddot{u} \) は疎行列ベクトル積）の当該方向の並進自由度成分のみを、当該層以上に属する節点について 集計し、符号を反転した値です。

\[ Q_i = -\sum_{j \geq i} \sum_{\text{dof} \in \text{階} j \text{ の並進自由度}} f_{\text{abs, dof}} \]

一貫質量行列では \( M \) が並進・回転自由度間の連成項を持ちますが、これは疎行列ベクトル積 \( M \cdot \ddot{u} \) の並進成分の中に既に正しく反映されているため、回転自由度そのものを 集計へ含める必要はありません（回転自由度の加速度に地動加速度 \( \ddot{x}_g \) を直接 加算することもしません）。どの階にも属さない節点（基礎レベルの自由節点等）の並進自由度は、 Z 座標が最も近い階（階の代表 Z＝所属節点の平均 Z）に計上します。最下層の \( Q_1 \) は 全並進自由度の慣性力の総和（ベースシア）に一致します。

層せん断力係数

全時刻を通じた層せん断力の絶対値最大値を、当該層以上の地震用重量の総和で除した値です。

\[ C_i = \dfrac{\max_t |Q_i(t)|}{\sum_{j \geq i} W_j} \]

\( W_j \) は階 \( j \) の地震用重量（5.1 地震力と同じ Story::seismic_weight）です。

階絶対加速度・階速度・階変位

階に属する節点の質量加重平均です（質量が全節点ゼロの階は単純平均）。加速度は絶対加速度（相対加速度＋地動加速度）、速度・変位は相対値です。

\[ \bar{a}_i = \dfrac{\sum_{\text{dof} \in \text{階} i} m_{\text{dof}} \left( \ddot{u}_{\text{dof}} + \ddot{x}_g \right)}{\sum_{\text{dof} \in \text{階} i} m_{\text{dof}}} \]

層応答の最大値

層せん断力・階絶対加速度・階速度・階変位のいずれも、全時刻を通じた絶対値最大値 （5.5.5 の間引きに関わらず全ステップで求めた値）を別途保持します。

実装：dynamic::timehistory（timehistory/recording.rs::ThRecorder）が算定します。X・Y 方向の両方を常に記録し、加振されていない方向の連成応答（ねじれ・偏心による応答等）も含みます。

5.5.5 詳細記録のフレーム間引き（record_every）

5.5.1〜5.5.2 いずれの経路でも、全節点変位・部材内力・層応答の時系列（5.5.4）は、 解析の全ステップのうち一定間隔 record_every ステップごとに間引いて保持します （最終ステップは間引きに関わらず必ず含めます）。地震波は数千〜数万ステップに及ぶため、 時系列の全量をそのまま保持するとメモリ・保存容量を圧迫することへの対処です。

record_every を指定しない場合は、記録フレーム数がおおむね 1000 になるよう次式で自動決定します（\( N \) は解析の全ステップ数）。

\[ \text{record_every} = \max\left(1, \left\lfloor N / 1000 \right\rfloor\right) \]

間引きの対象は時系列本体（各フレームの全節点変位・部材内力・層応答）のみです。 全ステップを通じた最大応答（節点変位の最大値、部材内力の包絡値、層せん断力係数 Ci、 層せん断力・階絶対加速度・階速度・階変位の絶対値最大）は間引かず、全ステップを走査して 求めます。

フレームごとの部材内力（時系列本体）は、各部材の評価断面のうち両端 2 点（材端、 最小ξ・最大ξ）のみを保持し、メモリを削減します（3D アニメーション・履歴表示は端部値の みを用いるため）。中間の評価断面（危険断面）を含む全断面の情報は、間引かない包絡値 （部材内力の最大値）側にのみ保持されます。

解析タブ「時刻歴応答」の「記録間引き」欄で record_every を指定できます （線形 Newmark-β・非線形のいずれの経路でも共通）。0 を指定すると上式による 自動決定になります。

実装：timehistory/recording.rs::auto_record_every が自動決定を、 呼び出し側（GUI・非線形の NonlinearThCfg::record_every）が record_every を 明示指定できます。