衛星構体の振動モード解析 — 打ち上げ環境とFEMによる予測

打ち上げの数分間で衛星にかかる加速度・振動・音響は、軌道上で衛星が一生のうちに経験するどんな擾乱よりも遥かに激しいものです。実際、宇宙開発の歴史を振り返ると「軌道に乗った瞬間に通信が途絶する」「電源系のリレーが脱落していた」「太陽電池パドルの固縛機構が壊れていた」といった事故の多くは、設計時の振動環境見積もりの甘さに端を発しています。フェアリング内の機器が「鳴き」、ロケットエンジンの低周波振動と衛星1次モードが共振してしまえば、設計加速度の何倍もの応答が局所的に発生してハーネスが断線します。

ではエンジニアは、打ち上げ環境を実機で再現することなく、どうやって「壊れない衛星」を設計するのでしょうか。鍵は モード解析（modal analysis） と FEM（有限要素法） にあります。衛星構体を質量と剛性の集合体として離散化し、その固有振動数と固有モードを事前に求めてしまえば、打ち上げロケットが規定する正弦・ランダム・衝撃環境に対する応答を「重ね合わせ」で予測できます。設計段階で「1次モードを横30Hz以上にする」「主構造の応答倍率を10倍以下に抑える」といった定量的な目標を立てられるようになるのです。

この技術は身近な応用先を多数持ちます。たとえば、CubeSat規格 (CDS) では構体1次モードの下限がデプロイヤと整合するよう規定されており、ベンダーは打ち上げ前の振動試験で固有振動数を確認します。Falcon 9 のユーザーマニュアルは衛星1次モードを「横方向 ≥ 25 Hz、縦方向 ≥ 35 Hz」と要求し、これを満たさない衛星は搭載できません。さらに、JWST のような大型観測衛星では、サブシステムの数百個の機器それぞれが振動応答スペクトルに対して個別認証されます。

本記事の内容

• 打ち上げ時に衛星構体が受ける振動環境（準静的、サイン、ランダム、音響、衝撃）の整理
• 多自由度系のモード解析の数学的導出 — 固有値問題、質量・剛性直交性、有効質量
• FEM による剛性・質量行列の組立と、Lanczos 法などの大規模固有値ソルバの考え方
• ロケット側 Coupled Load Analysis と衛星1次固有振動数の設計基準
• Python で「4自由度集中質量モデル」「Euler–Bernoulli 梁の連続体モード」「衝撃応答スペクトル(SRS)」を実装

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• ニュートンの運動の3法則を数式で理解する
• 単振動と減衰振動の理論と実装
• 連成振動と固有モードの解析
• 梁のたわみと弾性曲線方程式をわかりやすく解説
• 有限要素法（FEM）入門を基礎から徹底解説

直感: なぜ衛星は振動で壊れるか

衛星はロケットの先端にあるフェアリング内に固定され、打ち上げの最初の8分ほどで「自重の数倍の準静的加速度」「数十Hzから数千Hzの振動」「機体外壁を叩く140デシベル超の音響」という3種類の負荷を同時に受けます。仮に衛星の質量が500 kg で、ピーク時の加速度が 6 G ならば、構体には常時 30,000 N の慣性力が作用し続けることになります。これは小型自動車1台分を構体上面に静かに乗せ続けるのに近い負荷です。

しかし「壊れる」原因の多くは静的な大荷重ではなく、動的な共振応答 にあります。衛星構体のような連続体には固有振動数 $\omega_1, \omega_2, \dots$ が無数に存在し、外力の周波数がこれらに一致したとき、応答振幅は 減衰比 $\zeta$ に反比例する形 で増幅されます。典型的な衛星構体の減衰比は $\zeta \sim 0.01$（1%）程度で、共振時の応答倍率（Q値）は $Q = 1/(2\zeta) \sim 50$ にも達します。つまり入力 1 G の振動が共振帯では 50 G として伝わるのです。

ロケットの打ち上げ加速度プロファイルにはエンジン低周波（5–50 Hz）、推進系の脈動（数十–数百Hz）、空力騒音由来の広帯域成分（数百–数千 Hz）が混在しています。衛星の1次曲げモードがこの帯域のどこに乗るかが、設計の生死を分けます。1次モードを 十分高い周波数（30–40 Hz以上）に追い出す ことが、もっとも基本的な設計指針となります。

ここで自然な疑問が生まれます。「何種類の振動がどの帯域に存在するのか」を整理せずに固有振動数の話は始められません。次節で打ち上げ環境の振動分類を見ていきましょう。

打ち上げ環境の振動分類

ロケットメーカーが配布する ペイロードユーザーズマニュアル (PUG, Payload User’s Guide) には、衛星が経験する力学環境が階層的に整理されています。実機試験で再現可能な抽象化として、おおむね以下の5つに分けられます。

準静的加速度 (Quasi-Static, QSL)

リフトオフ時のスラスト立ち上がり、最大動圧 (Max-Q) 通過時のジャーク、段分離時の急減速など、比較的低周波（< 5 Hz）または直流的に作用する加速度です。Falcon 9 では縦方向に最大 +6.0 G / -2.0 G、横方向に ±2.0 G が QSL として規定されており、衛星構体はこれらを材料の降伏応力以下で支える必要があります。QSL は 静的解析（線形弾性 FEM） で検証され、安全率は通常 1.25（降伏）/ 1.4（破断）が用いられます。

サイン振動 (Sinusoidal Vibration, 5–100 Hz)

エンジン推力振動 (POGO)、機体軸振動、段間スロッシングなどに起因する 低周波の正弦振動 です。典型的なスペクトルは 5–10 Hz で 0.5 G、10–60 Hz で 1.0 G、60–100 Hz で減衰、といった矩形に近い包絡線で与えられます。サイン試験は基本的に 掃引（sweep）正弦波 で行い、5–100 Hz を毎オクターブ2–4分の速度で掃引しながら衛星基部加速度を制御します。この帯域は衛星の 1次・2次曲げモード（典型値 30–80 Hz）に直接ヒットするため、共振応答が最大の関心事になります。

ランダム振動 (Random Vibration, 20–2000 Hz)

エンジン乱流、空力騒音、機械系のフラッタなどがフェアリング内壁を励振し、広帯域のガウシアン的ランダム振動として機器搭載面に伝達されます。仕様は パワースペクトル密度 (PSD) で与えられ、Falcon 9 では 20 Hz で 0.005 $\mathrm{G^2/Hz}$、80–500 Hz で 0.04 $\mathrm{G^2/Hz}$ の plateau、2000 Hz で 0.005 $\mathrm{G^2/Hz}$ といった台形プロファイルが典型例です。RMS 加速度はおおむね 6–8 $G_\mathrm{rms}$ となり、3σ 値で 20 G を越えます。ランダム振動はサブシステム内部の電子部品（リレー、コンデンサ、PCB はんだ接合）の疲労破壊の主因です。

音響振動 (Acoustic Vibration, 140 dB)

リフトオフ時のエンジン排気音波と上昇中の空力騒音は、フェアリング内部に 130–145 dB の音場 を形成します。1/3オクターブバンドで指定され、全体音圧レベル (OASPL) で 140 dB を超えるのが典型です。音響は 大面積で軽量な構造（太陽電池パドル、アンテナ反射鏡、断熱多層膜 MLI）に対して、機械的ベース振動よりも遥かに激しい応答を引き起こします。設計上は 音響試験室（リバーブチャンバ） で実機を曝露します。

衝撃 (Pyro Shock, SRS)

段分離火工品、フェアリング分離、衛星分離（クランプバンド、ロックリング）といったパイロテクニックイベントは、ミリ秒オーダーで数千 G に達する 過渡的高周波衝撃 を発生させます。衝撃は単一の時間波形ではなく 衝撃応答スペクトル (Shock Response Spectrum, SRS) で仕様化されます。SRS は「各固有振動数 $f$ を持つ1自由度ばね–質量系に同じ時間波形を入力したとき、それぞれが示す最大応答」をプロットしたものです。Falcon 9 のクランプバンド分離 SRS は 100 Hz で 30 G、1 kHz で 1000 G、10 kHz で 4000 G といったプロファイルが典型です。

5つの環境のうち、サイン・ランダム・衝撃の3つはすべて 「外力の周波数に対する構造の応答」 を主題にしています。応答の予測には、まず衛星構体そのものの 固有振動数と固有モード が分かっていなければなりません。次節ではその数学的定式化を見ていきます。

モード解析の数学的定式化

自由振動の運動方程式

衛星構体を有限個の質量点と弾性ばね（あるいは弾性要素）の集合体として離散化し、自由度 $n$ のベクトル $\bm{u}(t) \in \mathbb{R}^n$ で各点の変位を表します。減衰を無視した自由振動の運動方程式は、ニュートン第二法則の多自由度版として次のように書けます。

$$ \bm{M} \ddot{\bm{u}}(t) + \bm{K} \bm{u}(t) = \bm{0} $$

ここで $\bm{M} \in \mathbb{R}^{n \times n}$ は 質量行列（正定値対称）、$\bm{K} \in \mathbb{R}^{n \times n}$ は 剛性行列（半正定値対称）です。1自由度の単振動 $m \ddot{u} + k u = 0$ の自然な拡張であり、各自由度が他の自由度とばねを通じて結合している様子を行列の非対角成分が表現しています。

固有値問題への帰着

「自由振動下で構造はどんな運動をするか」を知るには、全自由度が同じ周波数で同位相に振動する解（基準振動、normal mode）を仮定します。

$$ \bm{u}(t) = \bm{\phi}\, e^{j\omega t} $$

ここで $\bm{\phi} \in \mathbb{R}^n$ は時間に依存しないモード形状ベクトル、$\omega$ は角振動数です。これを運動方程式に代入すると、

$$ -\omega^2 \bm{M} \bm{\phi}\, e^{j\omega t} + \bm{K} \bm{\phi}\, e^{j\omega t} = \bm{0} $$

両辺を $e^{j\omega t} \neq 0$ で割って、

$$ \boxed{\,(\bm{K} – \omega^2 \bm{M})\, \bm{\phi} = \bm{0}\,} $$

を得ます。これは 一般化固有値問題 (generalized eigenvalue problem) で、$\bm{\phi} \neq \bm{0}$ となる非自明解を持つ条件は $\det(\bm{K} – \omega^2 \bm{M}) = 0$、すなわち $n$ 次の代数方程式です。この特性方程式の根として、$n$ 個の固有値 $\omega_i^2$（$i = 1, \dots, n$）と対応する固有ベクトル $\bm{\phi}_i$ が得られます。物理的な並びとして $\omega_1 \leq \omega_2 \leq \dots \leq \omega_n$ と昇順に並べ、$f_i = \omega_i / (2\pi)$ を $i$ 次固有振動数 [Hz] と呼びます。

質量直交性と剛性直交性

固有モード $\bm{\phi}_i, \bm{\phi}_j$（$i \neq j$）は、質量行列と剛性行列に関して 直交 します。これを示すには、二つの固有値方程式を並べます。

$$ \bm{K} \bm{\phi}_i = \omega_i^2 \bm{M} \bm{\phi}_i, \qquad \bm{K} \bm{\phi}_j = \omega_j^2 \bm{M} \bm{\phi}_j $$

第1式の両辺に左から $\bm{\phi}_j^\top$ をかけると $\bm{\phi}_j^\top \bm{K} \bm{\phi}_i = \omega_i^2 \bm{\phi}_j^\top \bm{M} \bm{\phi}_i$。同様に第2式の両辺に左から $\bm{\phi}_i^\top$ をかけて $\bm{\phi}_i^\top \bm{K} \bm{\phi}_j = \omega_j^2 \bm{\phi}_i^\top \bm{M} \bm{\phi}_j$ を得ます。$\bm{M}, \bm{K}$ が対称なので $\bm{\phi}_j^\top \bm{K} \bm{\phi}_i = \bm{\phi}_i^\top \bm{K} \bm{\phi}_j$ かつ $\bm{\phi}_j^\top \bm{M} \bm{\phi}_i = \bm{\phi}_i^\top \bm{M} \bm{\phi}_j$ ですから、両式を辺々引くと

$$ 0 = (\omega_i^2 – \omega_j^2)\, \bm{\phi}_i^\top \bm{M} \bm{\phi}_j $$

となります。固有値が異なる ($\omega_i \neq \omega_j$) という前提のもとで、

$$ \bm{\phi}_i^\top \bm{M} \bm{\phi}_j = 0 \quad (i \neq j) $$

すなわち 質量直交性 が成立します。これを上の式に代入すると同時に $\bm{\phi}_i^\top \bm{K} \bm{\phi}_j = 0$（剛性直交性）も従います。さらにモードを 質量正規化 (mass normalization) すれば、

$$ \bm{\phi}_i^\top \bm{M} \bm{\phi}_j = \delta_{ij}, \qquad \bm{\phi}_i^\top \bm{K} \bm{\phi}_j = \omega_i^2 \delta_{ij} $$

という美しい形にまとまります。ここで $\delta_{ij}$ はクロネッカーのデルタです。

モード重ね合わせ

直交性が成立すると、任意の運動 $\bm{u}(t)$ をモードベクトルの線形結合で展開できます。

$$ \bm{u}(t) = \sum_{i=1}^{n} q_i(t)\, \bm{\phi}_i $$

ここで $q_i(t)$ は モード座標 (modal coordinate) あるいは 基準座標 と呼ばれるスカラー量です。これを運動方程式（外力 $\bm{F}(t)$ と Rayleigh 減衰 $\bm{C}$ も含めて）

$$ \bm{M}\ddot{\bm{u}} + \bm{C}\dot{\bm{u}} + \bm{K}\bm{u} = \bm{F}(t) $$

に代入し、左から $\bm{\phi}_i^\top$ をかけて直交性を使うと、$n$ 本の 独立な1自由度方程式 に分離されます。

$$ \ddot{q}_i + 2\zeta_i \omega_i \dot{q}_i + \omega_i^2 q_i = \bm{\phi}_i^\top \bm{F}(t) $$

ここで $\zeta_i$ は $i$ 次モードの減衰比です。多自由度連成系の解析が、$n$ 個の独立な単振動子の重ね合わせに帰着するという驚くべき結果が得られました。サイン応答もランダム応答も SRS も、すべてこの「モードごとの単振動子応答」の重ね合わせとして計算できるのです。

有効質量

設計実務で重要になるのが 有効質量 (effective modal mass) $M_{\mathrm{eff},i}$ です。基部加速度 $\ddot{u}_b$ が一様に作用するとき、$i$ 次モードに参加する質量分が $M_{\mathrm{eff},i}$ で、

$$ M_{\mathrm{eff},i} = \frac{(\bm{\phi}_i^\top \bm{M} \bm{r})^2}{\bm{\phi}_i^\top \bm{M} \bm{\phi}_i} $$

と定義されます。ここで $\bm{r}$ は 影響ベクトル で、基部の単位変位がもたらす各自由度の変位を表します（並進では全要素1のベクトル）。総質量保存則として $\sum_i M_{\mathrm{eff},i} = M_\mathrm{tot}$ が成り立ちます。実務では「累積有効質量が全質量の 90% 以上に達するまで」モードを抽出することが標準的な判定基準です。

ここまで、自由振動の運動方程式から固有値問題、モード重ね合わせ、有効質量までを導きました。次に、実際の衛星構体のような複雑な形状で $\bm{M}, \bm{K}$ をどう構築し、固有値問題をどう解くかを見ていきます。

固有値問題の解法 (FEM)

要素分割と行列組立

実際の衛星構体は曲がった円筒、ハニカムサンドイッチパネル、ハーネスや機器が点在する不規則な3次元形状をしています。これを 有限要素 (finite element) に分割し、各要素の弾性・慣性を行列で表現してから、全体の $\bm{M}, \bm{K}$ を組み立てるのが FEM の基本戦略です。

衛星構体で典型的に使われる要素タイプは:

• 梁要素 (beam element) — トラスアセンブリ、機器搭載ブラケット、ロックリング
• シェル要素 (shell element) — ハニカムパネル、円筒タンク壁、太陽電池基板
• 固体要素 (solid element) — 接手部、機器筐体、推進系バルブ取付部
• 集中質量要素 (point mass) — リアクションホイール、バッテリ、推進剤

各要素の質量行列 $\bm{M}^e$ と剛性行列 $\bm{K}^e$ は要素の形状関数と弾性係数から導かれます。たとえば長さ $L_e$、断面積 $A_e$、密度 $\rho$、ヤング率 $E$ の1次元バー要素では、整合質量行列と剛性行列がそれぞれ

$$ \bm{M}^e = \frac{\rho A_e L_e}{6}\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}, \qquad \bm{K}^e = \frac{E A_e}{L_e}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} $$

となります（詳しい導出は 有限要素法（FEM）入門 を参照）。全体の $\bm{M}, \bm{K}$ は各要素の節点を全体節点番号にマップしながら、対応する位置に要素行列を加算（assemble）して構築します。

境界条件の処理

衛星はロケット側の PAF (Payload Attach Fitting) にクランプバンド等で固定されます。モデル上では「PAF 部の節点を完全固定」と扱い、対応する自由度を $\bm{M}, \bm{K}$ から削除（あるいはペナルティ法で大剛性化）します。これにより固有値問題が 正定値 となり、$\omega_1 > 0$ が得られます（自由境界ではゼロ固有値 = 剛体モード6つが現れます）。

大規模固有値ソルバ

衛星 FEM の自由度数は数十万から数百万に達するため、$n \times n$ 行列をフルに対角化することは現実的ではありません。実務的には 下から数十モード だけを高効率に抽出するアルゴリズムが用いられます。

• 逆反復法 (inverse iteration) — 最小固有値に対応するベクトルへの収束。$\bm{K}^{-1} \bm{M} \bm{\phi}$ を反復。
• シフト付き逆反復 — 指定周波数近傍のモードを抽出。
• Lanczos 法 — 三重対角行列に縮約してから内部固有値問題を解く。数百モードを効率的に抽出可能。
• 部分空間反復 (subspace iteration) — 複数の Ritz ベクトルを同時に直交化しながら反復。

商用 FEM ソルバ（NASTRAN、ABAQUS、ANSYS）ではほぼ全てが内部で Lanczos 法ベースのソルバを採用しています。実用上、解析者がこれらの数値手法を直接実装することは稀ですが、「なぜ FEM の固有値解析は『小さな固有値から順に』取れるのか」を理解する上で、逆反復のイメージは重要です。

解析精度と要素サイズ

固有振動数の収束精度は要素サイズ $h$ に依存し、典型的には $O(h^2)$ の収束次数を持ちます。目安として「対象モードの半波長あたり 4–6 要素」 が経験則です。たとえば 100 Hz の曲げモードの半波長が 1 m なら、200 mm 程度のメッシュサイズが必要です。これより粗いと固有振動数を 過大評価 する（剛すぎる近似になる）傾向があります。

ここまでで「衛星構体の固有値を FEM でどう求めるか」が見えました。次は、ロケット側がどんな固有振動数を要求しているか、その設計基準を確認しましょう。

設計基準と1次固有振動数

Coupled Load Analysis (CLA)

ロケット会社は、自社の機体ダイナミクスモデルと顧客衛星のダイナミクスモデル（縮約モデル）を結合し、各打ち上げイベント（リフトオフ、Max-Q、段分離、エンジンカットオフ）に対する 時間応答解析 を行います。これを Coupled Load Analysis (CLA) と呼びます。CLA の入力として、衛星側は Craig-Bampton 縮約モデル や 質量・剛性行列のサブセット をロケット側に提出する必要があります。

CLA から得られる出力は、衛星基部での加速度時刻歴・パワースペクトル・等価QSL です。これらが PUG の規定値（最大期待飛行レベル, MPE）を上回らないことが、打ち上げ許可の条件となります。

1次固有振動数の下限要求

CLA を成立させるためには、衛星と機体のモードが分離している ことが大前提です。具体的には、衛星の1次曲げモード（横方向）と1次縦モード（軸方向）が、ロケットの主要な低周波モードや POGO 周波数から十分離れている必要があります。これを担保するのが、

といった 下限値 です（表は近年公開値の代表例で、実際の契約条件は個別に決まります）。CubeSat のような小型衛星はデプロイヤ (P-POD, ISIPOD) に格納されるため、デプロイヤ自体の固有モードが要求になり、しばしば 90 Hz 以上といった高い値が課されます。

「30/40 Hz ルール」の本質

横30 / 縦40 という数値は経験則のように見えますが、その背後には 「衛星モードがロケットのエンジン推力振動帯（10–25 Hz）の上に出る」 という物理的要請があります。仮に衛星1次モードが 20 Hz にあり、ロケットエンジンが 18 Hz で振動していると、共振応答倍率 Q ~ 50 が掛かって衛星基部に 50 G の加速度が立ちます。これを 30 Hz 以上に追いやれば、エンジン振動帯から離れて Q 倍は劇的に下がり、設計加速度内に収まります。

設計フローへの組み込み

このため、衛星設計の初期段階（コンセプト設計）から、1次固有振動数を 40–60 Hz 以上（マージンを含む）に設定するのが標準的な習慣です。具体的な手法としては:

• 主構造を 円筒（モノコック）構造 とし、座屈と曲げ剛性を同時に確保する
• ハニカムサンドイッチパネル で軽量かつ高剛性を実現
• 重量機器（リアクションホイール、推進系タンク）を 基部近く に配置して有効モード質量を上げる
• 太陽電池パドルやアンテナなど 柔らかい付属物は固縛機構 で打ち上げ時に剛体化

衛星のシステムエンジニアと構造エンジニアは、コンセプト段階からこれらの要求を相互に調整します。

設計基準を満たしているかは、最終的には FEM で1次モードを計算して確認します。続いて、Python で最も単純な集中質量モデルから FEM まで段階的にモード解析を実装していきます。

Python実装 — 集中質量〜FEMモード解析

4自由度集中質量モデル

最初に、衛星構体を「4段の集中質量がばねで連結された塔」として極度に単純化したモデルでモード解析の流れを掴みます。下から PAF、構体下段、構体上段、機器プラットフォームの順に質量 $m_1, \dots, m_4$ がばね定数 $k_1, \dots, k_4$ で結合されていると考えます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

# 4自由度集中質量モデル(下から順)
m = np.array([50.0, 80.0, 100.0, 200.0])      # [kg] 各層質量
k = np.array([3.0e7, 2.0e7, 1.5e7, 1.0e7])    # [N/m] 各層ばね定数

# 質量行列（対角）
M = np.diag(m)

# 剛性行列の組立(隣接層間ばねが両端に対称に効く)
n = len(m)
K = np.zeros((n, n))
for i in range(n):
    K[i, i] += k[i]
    if i + 1 < n:
        K[i, i]   += k[i + 1]
        K[i+1, i+1] += 0.0       # 後で足す
        K[i, i+1] -= k[i + 1]
        K[i+1, i] -= k[i + 1]
# 最上層には外向きばねがないので k[3] を1回しか足していない
# 修正: 上端を自由端として処理
K[n-1, n-1] = k[n-1]
K[n-1, n-2] = -k[n-1]
K[n-2, n-1] = -k[n-1]

print("Mass matrix M =\n", M)
print("Stiffness matrix K =\n", K)

# 一般化固有値問題 K phi = w^2 M phi
eigvals, eigvecs = eigh(K, M)
omega = np.sqrt(eigvals)        # 角振動数 [rad/s]
freq  = omega / (2 * np.pi)     # 周波数 [Hz]

print("\n固有振動数 [Hz]:", np.round(freq, 2))

このコードでは scipy.linalg.eigh を使って一般化固有値問題 $\bm{K}\bm{\phi} = \omega^2 \bm{M}\bm{\phi}$ を直接解いています。eigh は対称正定値ペアに対する高速かつ数値安定なルーチンで、内部的には Cholesky 分解後に標準固有値問題に帰着させる方式を採ります。出力される eigvals は $\omega^2$ の昇順、eigvecs の各列が対応する質量正規化された固有モード $\bm{\phi}_i$ です。実行すると、最低次モードはおおむね 18 Hz、最高次は 200 Hz 程度の値になります。質量の重い最上段（200 kg）を低周波で揺らす1次モードと、軽い最下段が高周波で振動する4次モードという質量分布と振動数の対応が直感的に確認できます。

続いて、得られたモード形状を可視化します。

import numpy as np
import matplotlib.pyplot as plt

heights = np.array([0.5, 1.0, 1.5, 2.0])    # 各質量点の高さ [m]
# 基部 (0 m) を加えてプロット
z = np.concatenate([[0.0], heights])

fig, axes = plt.subplots(1, n, figsize=(12, 4.5), sharey=True)
for i in range(n):
    phi = np.concatenate([[0.0], eigvecs[:, i]])
    phi = phi / np.max(np.abs(phi))           # ピーク振幅で正規化
    axes[i].plot(phi, z, 'o-', lw=2)
    axes[i].axvline(0, color='gray', alpha=0.5, ls='--')
    axes[i].set_xlim(-1.2, 1.2)
    axes[i].set_title(f'Mode {i+1}\n{freq[i]:.1f} Hz')
    axes[i].set_xlabel('Normalized displacement')
    axes[i].grid(alpha=0.3)
axes[0].set_ylabel('Height [m]')
plt.suptitle('Modal shapes of 4-DOF satellite model')
plt.tight_layout()
plt.savefig('modal_shapes_4dof.png', dpi=150, bbox_inches='tight')
plt.show()

得られた4本のグラフから、モード形状の物理的な意味が明瞭に読み取れます。1次モードは全質量が同位相で大きく揺れる「片持ち梁の1次曲げ」型で、節（ゼロ交差点）が基部にしかありません。2次モードは中段に節が1つ現れ、上下が逆位相に揺れます。3次・4次となるほど節の数が増えていき、空間的に短波長の振動になります。これは弦の倍音や Euler–Bernoulli 梁の高次モードとまったく同じ性質で、自由度を増やせば連続体のモードに収束します。

有効モード質量の計算

設計判定で必要となる有効モード質量を計算します。

import numpy as np

# 影響ベクトル(基部一様加速度に対して、全質量が同じ単位変位を取る)
r = np.ones(n)

# 質量正規化(eighの出力は既に質量正規化されている)
# M_eff_i = (phi_i^T M r)^2 / (phi_i^T M phi_i)
M_eff = np.zeros(n)
for i in range(n):
    phi = eigvecs[:, i]
    num = (phi @ M @ r) ** 2
    den = phi @ M @ phi
    M_eff[i] = num / den

ratio = M_eff / m.sum() * 100
cum   = np.cumsum(ratio)

print("\n--- Modal mass participation ---")
for i in range(n):
    print(f"Mode {i+1}: f={freq[i]:5.1f} Hz, "
          f"M_eff={M_eff[i]:6.1f} kg ({ratio[i]:5.1f}%), "
          f"cumulative {cum[i]:5.1f}%")
print(f"Total mass = {m.sum():.1f} kg, sum of M_eff = {M_eff.sum():.1f} kg")

このコードを実行すると、典型的には1次モードに全質量の70–85%が、2次モードにさらに10–20%が集まり、累積で2次までに95%以上に達することがわかります。これが意味するのは、「外部加速度に対する応答の大半は1–2次モードだけで決まる」という重要な事実です。実機 FEM でも、有効質量の累積が90%を超えるまでのモード数を抽出すれば、CLA に十分な精度が得られるとされています。逆に「高次モードの固有振動数は分かっているが有効質量が極端に小さい」モードは、設計影響が小さいので無視できます。

Euler–Bernoulli梁による連続体モード

集中質量モデルは直感的ですが、衛星構体の実態は「連続体としての梁・シェル」です。連続体の解析的解との比較は、FEM 結果の妥当性検証 (validation) で重要になります。長さ $L$、曲げ剛性 $EI$、線密度 $\rho A$ の片持ち梁 (cantilever) の自由振動方程式

$$ \rho A \frac{\partial^2 w}{\partial t^2} + EI \frac{\partial^4 w}{\partial x^4} = 0 $$

を $w(x, t) = \phi(x) e^{j\omega t}$ で変数分離すると、$\phi(x)$ は4階常微分方程式 $\phi”” = \beta^4 \phi$ を満たします（ここで $\beta^4 = \rho A \omega^2 / EI$）。境界条件「基部固定 ($\phi=\phi’=0$)、自由端自由 ($\phi”=\phi”’=0$)」を課すと、特性方程式

$$ \cos(\beta L) \cosh(\beta L) + 1 = 0 $$

の根 $\beta_i L$ から $i$ 次固有振動数 $\omega_i = (\beta_i L)^2 \sqrt{EI/(\rho A L^4)}$ が決まります。最初の数根の数値は $\beta_1 L = 1.8751$, $\beta_2 L = 4.6941$, $\beta_3 L = 7.8548$ です。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq

# 衛星アダプタ円筒を模擬した諸元
L  = 1.5                      # [m] 長さ
E  = 70e9                     # [Pa] アルミ
rho = 2700.0                  # [kg/m^3]
D_out, D_in = 0.40, 0.396     # [m] 外径/内径(薄肉円筒)
A  = np.pi/4 * (D_out**2 - D_in**2)
I  = np.pi/64 * (D_out**4 - D_in**4)
EI  = E * I
rhoA = rho * A
print(f"Mass per length = {rhoA:.2f} kg/m, total mass = {rhoA*L:.1f} kg")
print(f"EI = {EI:.3e} N m^2")

# 特性方程式 cos(bL)cosh(bL) + 1 = 0 の根を数値的に求める
def char_eq(bL):
    # 大きい bL で cosh が発散するのを避けるため正規化
    return np.cos(bL) + 1.0 / np.cosh(bL)

# 最初の5モードの bL を求める(根の近傍ブラケットを与える)
bL_roots = []
brackets = [(1.0, 3.0), (4.0, 5.5), (7.0, 9.0), (10.0, 11.5), (13.0, 14.5)]
for lo, hi in brackets:
    bL_roots.append(brentq(char_eq, lo, hi))
bL_roots = np.array(bL_roots)
print("\nbeta_i L =", np.round(bL_roots, 4))

# 角振動数と周波数
omega_analytic = bL_roots**2 * np.sqrt(EI / (rhoA * L**4))
freq_analytic  = omega_analytic / (2*np.pi)
print("固有振動数 [Hz] (analytic) =", np.round(freq_analytic, 2))

ここで使ったテクニックは、$\cos(\beta L)\cosh(\beta L) + 1 = 0$ を $\cos(\beta L) + 1/\cosh(\beta L) = 0$ と書き換える正規化です。$\beta L$ が10を超えると $\cosh$ は数万倍に発散して数値オーバーフローを引き起こすので、両辺を $\cosh(\beta L)$ で割っておくと安定します。実行すると、長さ1.5 m、外径400 mm のアルミ薄肉円筒の1次固有振動数はおおむね 100 Hz 前後に出ます。CubeSat ほど小型なら数百 Hz、大型衛星なら数十 Hz が典型値です。

FEM 梁モデルとの比較

同じ円筒を FEM の Euler–Bernoulli 梁要素でモード解析し、解析解との一致を確認します。1要素あたりに節点を2つ持ち、各節点に「たわみ $w$」と「回転 $\theta$」の2自由度を割り当てる Hermite 3次補間 が標準です。要素の長さを $h$ とすると、整合質量・剛性行列は次のように与えられます。

$$ \bm{K}^e = \frac{EI}{h^3}\begin{bmatrix} 12 & 6h & -12 & 6h \\ 6h & 4h^2 & -6h & 2h^2 \\ -12 & -6h & 12 & -6h \\ 6h & 2h^2 & -6h & 4h^2 \end{bmatrix} $$

$$ \bm{M}^e = \frac{\rho A h}{420}\begin{bmatrix} 156 & 22h & 54 & -13h \\ 22h & 4h^2 & 13h & -3h^2 \\ 54 & 13h & 156 & -22h \\ -13h & -3h^2 & -22h & 4h^2 \end{bmatrix} $$

この要素行列を組み立てて固有値を計算します。

import numpy as np
from scipy.linalg import eigh
import matplotlib.pyplot as plt

def beam_modal_fem(L, EI, rhoA, n_elem=20):
    """片持ち梁のFEMモード解析(Euler-Bernoulli, Hermite要素)"""
    h = L / n_elem
    n_node = n_elem + 1
    ndof = 2 * n_node           # 各節点に w, theta

    # 要素行列
    Ke_base = EI / h**3 * np.array([
        [ 12,    6*h, -12,    6*h],
        [ 6*h, 4*h*h, -6*h, 2*h*h],
        [-12,   -6*h,  12,   -6*h],
        [ 6*h, 2*h*h, -6*h, 4*h*h]])
    Me_base = rhoA * h / 420 * np.array([
        [156,   22*h,   54,  -13*h],
        [22*h, 4*h*h,  13*h, -3*h*h],
        [54,   13*h,  156,  -22*h],
        [-13*h,-3*h*h,-22*h, 4*h*h]])

    K = np.zeros((ndof, ndof))
    M = np.zeros((ndof, ndof))
    for e in range(n_elem):
        idx = [2*e, 2*e+1, 2*e+2, 2*e+3]
        for i in range(4):
            for j in range(4):
                K[idx[i], idx[j]] += Ke_base[i, j]
                M[idx[i], idx[j]] += Me_base[i, j]

    # 境界条件: 基部固定(節点0の w, theta = 0)
    free = list(range(2, ndof))
    Kf = K[np.ix_(free, free)]
    Mf = M[np.ix_(free, free)]

    eigvals, eigvecs = eigh(Kf, Mf)
    omega = np.sqrt(eigvals)
    return omega, eigvecs, free, ndof, h

omega_fem, modes_fem, free_dof, ndof, h = beam_modal_fem(L, EI, rhoA, n_elem=40)
freq_fem = omega_fem / (2*np.pi)

print("\n--- FEM vs analytical (first 5 modes) ---")
print(f"{'mode':<6}{'FEM [Hz]':<14}{'Analytical [Hz]':<18}{'error [%]'}")
for i in range(5):
    err = (freq_fem[i] - freq_analytic[i]) / freq_analytic[i] * 100
    print(f"{i+1:<6}{freq_fem[i]:<14.3f}{freq_analytic[i]:<18.3f}{err:+.3f}")

40要素 FEM と解析解の差は1次モードで 0.01% 未満、5次モードでも 0.5% 程度に収まります。これは Hermite 3次要素の収束次数が $O(h^4)$ と非常に高く、対象モードあたり数要素あれば十分な精度が出ることを示しています。逆に、要素数を 5 まで減らすと高次モードの誤差は数パーセント以上に拡大します。実機 FEM でも「対象モードの半波長あたり 4–6 要素」を確保することが重要です。

共振応答倍率の可視化

固有振動数だけでなく、外力周波数 $\Omega$ に対する応答倍率 $H(\Omega)$ を可視化することで、「ロケットのサイン振動帯がモードに刺さるとどうなるか」が直感的に理解できます。1自由度系の周波数応答関数は

$$ |H(\Omega)| = \frac{1}{\sqrt{(1 – r^2)^2 + (2\zeta r)^2}}, \quad r = \Omega / \omega_n $$

で与えられます。複数モードの応答を合成して描いてみましょう。

import numpy as np
import matplotlib.pyplot as plt

freqs_natural = freq_fem[:5]              # 主要5モード
zetas = 0.01 * np.ones(5)                 # 1%減衰

f_drive = np.logspace(0, 3, 2000)         # 1 Hz - 1 kHz
H_total = np.zeros_like(f_drive)

plt.figure(figsize=(10, 6))
for fi, zi in zip(freqs_natural, zetas):
    r  = f_drive / fi
    Hi = 1.0 / np.sqrt((1 - r**2)**2 + (2 * zi * r)**2)
    plt.loglog(f_drive, Hi, alpha=0.4, lw=1.0,
               label=f'mode {fi:.0f} Hz')
    H_total = np.maximum(H_total, Hi)    # 包絡

plt.loglog(f_drive, H_total, 'k-', lw=2.0, label='envelope')
plt.axhspan(1, 100, color='lightyellow', alpha=0.3)
plt.axvline(30, color='red', ls='--', alpha=0.5, label='30 Hz spec')
plt.xlabel('Driving frequency [Hz]')
plt.ylabel('Amplification |H|')
plt.title('Resonance amplification of cantilever beam (zeta=1%)')
plt.legend(loc='lower left', fontsize=9)
plt.grid(True, which='both', alpha=0.3)
plt.ylim(0.01, 200)
plt.tight_layout()
plt.savefig('resonance_amplification.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、共振応答の支配的な特徴が読み取れます。各モードの固有振動数で応答倍率は約 $1/(2\zeta) = 50$ 倍にピークし、その左右では急峻に低下します。1次モード（赤線右側のピーク）が「30 Hz スペック」より上にあれば、ロケットのサイン振動帯（5–30 Hz）では応答倍率が1以下に抑えられ、入力加速度がそのまま伝達されるだけになります。逆に、もし1次モードが 20 Hz にあれば、20 Hz の入力に対して 50 倍の応答が立ち、構体は容易に破壊します。これが「1次モードを高周波側に追い出す」設計指針の物理的根拠です。

ここまでで、サイン振動とランダム振動に対するモード応答の予測ができました。最後に、過渡的な衝撃に対する応答評価の標準ツールである SRS を実装します。

衝撃応答スペクトル(SRS)

SRSの定義

衝撃応答スペクトル (Shock Response Spectrum, SRS) は、衝撃時間波形 $a_b(t)$ をある周波数 $f_n$ を持つ1自由度ばね–質量系に 基部加速度 として入力したときの、質量の 絶対加速度の最大値 $|a_\mathrm{max}|$ を、周波数を変えながらプロットしたものです。

$$ \mathrm{SRS}(f_n) = \max_t |a_\mathrm{response}(t; f_n, \zeta)| $$

1自由度系の基部励起時の運動方程式は、相対変位 $z = u – u_b$ について

$$ \ddot{z} + 2\zeta \omega_n \dot{z} + \omega_n^2 z = -\ddot{u}_b $$

となり、これを数値積分すれば応答が得られます。SRS は「衝撃の中身（時間波形）」を抽象化し、「あらゆる固有振動数の構造がどれだけ揺らされるか」 を一望できる強力な可視化です。クランプバンド分離や火工品の衝撃仕様は、ほぼ例外なく SRS で与えられます。

Smallwoodのデジタルフィルタ法

SRS の標準的な計算法は Smallwood (1981) のデジタルフィルタ法 です。1自由度系の伝達関数を双線形変換 (Tustin 変換) ではなく インパルス不変変換 で離散化し、入力サンプリング周波数 $f_s$ に対して2タップの IIR フィルタとして実装します。連続系の伝達関数

$$ H(s) = \frac{-1}{s^2 + 2\zeta \omega_n s + \omega_n^2} $$

（出力は相対加速度 $\ddot{z}$）に対し、減衰固有角振動数 $\omega_d = \omega_n \sqrt{1-\zeta^2}$ を用いた離散化係数は次のようになります（詳細導出は Smallwood の原論文を参照）。

$$ b_0 = 1 – \frac{\omega_n}{\omega_d} e^{-\zeta\omega_n T} \sin(\omega_d T), \quad b_1 = -2 b_0 \cos(\omega_d T) e^{-\zeta\omega_n T},\, \dots $$

これを実装します。

import numpy as np

def srs_smallwood(accel_input, dt, freqs, zeta=0.05):
    """Smallwood digital filter method による SRS 計算
    accel_input: 入力加速度時系列(配列), dt: サンプリング周期 [s],
    freqs: SRSを評価する周波数配列 [Hz], zeta: 減衰比"""
    n = len(accel_input)
    srs_pos = np.zeros(len(freqs))    # 正側ピーク
    srs_neg = np.zeros(len(freqs))    # 負側ピーク

    for k, fn in enumerate(freqs):
        wn = 2 * np.pi * fn
        wd = wn * np.sqrt(1 - zeta**2)
        E  = np.exp(-zeta * wn * dt)
        K_ = wd * dt
        C, S = np.cos(K_), np.sin(K_)

        # Smallwood係数(絶対加速度SRS)
        b0 = 1 - E * S / K_
        b1 = 2 * (E * S / K_ - E * C)
        b2 = E**2 - E * S / K_
        a1 = -2 * E * C
        a2 = E**2

        # IIRフィルタを直接実装
        y = np.zeros(n)
        x = accel_input
        for i in range(2, n):
            y[i] = (b0*x[i] + b1*x[i-1] + b2*x[i-2]
                    - a1*y[i-1] - a2*y[i-2])
        srs_pos[k] = np.max(y)
        srs_neg[k] = -np.min(y)
    # maximax SRS(両側ピークの大きい方)
    return np.maximum(srs_pos, srs_neg)

この関数は、入力加速度時系列を周波数 $f_n$ の1自由度共振系に通したときの出力の絶対値最大を、各 $f_n$ について計算します。係数の式中の K_ = wd * dt は離散化周期と固有周期の比に対応し、これが小さいほど（サンプリングが密なほど）連続系の応答に近づきます。実用上は、SRS を評価したい最高周波数の 10倍以上のサンプリング周波数 が必要です。

簡易ハーフサイン衝撃の SRS

実際の火工衝撃波形は複雑な減衰振動ですが、まずは教科書的な ハーフサイン衝撃（持続時間 $\tau$、ピーク加速度 $A_0$）で SRS の振る舞いを確認します。

import numpy as np
import matplotlib.pyplot as plt

# 入力衝撃: half-sine pulse
fs = 200_000.0                     # サンプリング 200 kHz
dt = 1.0 / fs
tau = 0.001                        # 衝撃持続時間 1 ms
A0  = 1000.0                       # ピーク加速度 1000 G
t_total = 0.05                     # 50 ms 観測
t = np.arange(0, t_total, dt)
a_in = np.zeros_like(t)
mask = t < tau
a_in[mask] = A0 * np.sin(np.pi * t[mask] / tau)

# SRSを評価する周波数(対数等間隔, 10 Hz - 10 kHz)
freqs = np.logspace(1, 4, 60)
srs = srs_smallwood(a_in, dt, freqs, zeta=0.05)

fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
ax[0].plot(t*1000, a_in)
ax[0].set_xlim(0, 5)
ax[0].set_xlabel('Time [ms]')
ax[0].set_ylabel('Input acceleration [G]')
ax[0].set_title(f'Half-sine pulse: A0={A0} G, tau={tau*1000:.1f} ms')
ax[0].grid(alpha=0.3)

ax[1].loglog(freqs, srs, 'o-', lw=1.5)
ax[1].axhline(A0, color='gray', ls=':', label=f'Input peak = {A0} G')
ax[1].axvline(1/tau, color='red', ls='--', alpha=0.5,
              label=f'1/tau = {1/tau:.0f} Hz')
ax[1].set_xlabel('Natural frequency [Hz]')
ax[1].set_ylabel('SRS [G]')
ax[1].set_title('Shock Response Spectrum (zeta=5%)')
ax[1].legend()
ax[1].grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('srs_halfsine.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、SRS の特徴的な3つの帯域が明瞭に読み取れます。低周波領域（$f_n \ll 1/\tau$）では、衝撃の継続時間に対して系の固有周期が長すぎるため応答は小さく、入力速度変化に比例して $f_n$ に対して直線的に増加します。中間領域（$f_n \sim 1/\tau$）では、入力周波数成分と系の固有振動数が共鳴し、SRS は入力ピーク $A_0$ の約1.7倍程度まで増幅します。高周波領域（$f_n \gg 1/\tau$）では系が「瞬時に追従」して入力ピーク値に漸近します。火工衝撃の SRS 仕様が「10 Hz から 10 kHz にかけて低周波で 30 G、高周波で 1000 G 以上」のような 低周波小→高周波大 のスロープを持つのは、衝撃のエネルギーが短時間に集中している（高周波成分が大きい）ことの直接の現れです。

モード解析と SRS の結合 — 機器搭載部の応答

最後に、これまでに求めた構体の固有モード情報と SRS を結合して、「ロケット側から指定された SRS に対し、機器搭載点はどれだけの加速度を受けるか」を見積もります。簡易には「機器搭載点の応答 ≒ 構体1次固有振動数における SRS 値 × 機器搭載点のモード係数」と評価できます。

import numpy as np

# 構体1次固有振動数 (FEM結果)
f1 = freq_fem[0]
# 機器搭載点 (先端) のモード変位
phi_tip = modes_fem[-2, 0]    # 先端ノードのたわみ成分
# モード参加係数
gamma1 = np.sum(modes_fem[:, 0])    # 簡易近似
# 1次モードの SRS 値
srs_at_f1 = np.interp(f1, freqs, srs)
print(f"構体1次固有振動数: {f1:.1f} Hz")
print(f"SRS @ {f1:.1f} Hz : {srs_at_f1:.1f} G")
print(f"先端搭載機器の概算ピーク加速度: {srs_at_f1 * abs(phi_tip / modes_fem[0, 0]) :.1f} G")

この簡易計算は厳密ではありませんが、「SRS 上の構体共振点の値が、そのまま構体応答のオーダーを決める」という重要な感覚を与えてくれます。実機設計では、FEM で全モードを求めた後、各モードに対して SRS から励起レベルを読み取り、モード重ね合わせ（あるいは絶対値和、SRSS 和）で機器搭載点の応答を予測する手順を取ります。Falcon 9 の SRS が 100 Hz で 30 G しかなくても、衛星1次モードが 100 Hz にあり Q = 50 ならば、機器搭載部は 30 × 50 = 1500 G 級の応答を受け得るのです。

まとめ

本記事では、衛星構体の振動環境とモード解析、FEM による固有値解析、そして衝撃応答スペクトルまでを体系的に解説しました。

• 打ち上げ環境の5分類 — 準静的 (< 5 Hz)、サイン (5–100 Hz)、ランダム (20–2000 Hz)、音響 (140 dB)、衝撃 (SRS)。これらすべてが「外力周波数 × 構造応答」の枠組みで扱える。
• モード解析の数学的基盤 — 自由振動方程式 $\bm{M}\ddot{\bm{u}} + \bm{K}\bm{u} = \bm{0}$ から一般化固有値問題 $(\bm{K} – \omega^2 \bm{M})\bm{\phi} = \bm{0}$ が導かれ、質量・剛性直交性によりモード座標で系が独立な1自由度方程式に分解される。
• 有効モード質量 — $M_{\mathrm{eff},i} = (\bm{\phi}_i^\top \bm{M} \bm{r})^2 / (\bm{\phi}_i^\top \bm{M} \bm{\phi}_i)$ で定義され、累積90%までのモードが応答の大半を支配する。
• FEM 固有値解析 — 梁・シェル・固体要素で行列を組立て、Lanczos 法等で低次モードを抽出。要素サイズは対象モード半波長あたり 4–6 要素が目安。
• 設計基準 — Falcon 9 は横 25 / 縦 35 Hz、H3 は横 30 / 縦 40 Hz が衛星1次モードの下限。これによりロケットエンジン振動帯との分離を担保。
• 衝撃応答スペクトル — Smallwood のデジタルフィルタ法で計算され、火工衝撃の周波数依存応答を一望できる。共振が刺さると入力 SRS の Q 倍が機器搭載部に立つ。

これらの手法はすべて、設計の初期段階から振動環境を見抜き、「壊れる前に直す」ための工学的言語です。CubeSat から大型観測衛星まで、構造解析担当者が日常的に使う FEM ソルバの裏側で起きていることが、本記事の方程式そのものです。

次のステップとして、以下の記事も参考にしてください。

• 単振動と減衰振動の理論と実装
• 連成振動と固有モードの解析
• 有限要素法（FEM）入門を基礎から徹底解説
• 梁のたわみと弾性曲線方程式をわかりやすく解説
• オイラー座屈の理論をわかりやすく解説
• 曲げ応力と断面二次モーメントの関係