立体骨組構造の固有振動解析結果の検算
Check of Free Vibration Analysis of 3D Framed Structure

SGST 正会員 瀧上工業（株）

安 藤 浩 吉SGST 正会員 瀧上工業（株）技術部 工博 織 田 博 孝
名城大学理工学部教授 工博 中 川 建 治

まえがき
昨今の骨組構造の固有振動解析にはほとんどにコンピュータが利用され、入力データを用意すれば経験の少ない技術者でも容易に結果を得ることができる。このため、信頼できるプログラムを使用するにもかかわらず、使用法の誤り、入力データの意味の取違い、データ自身の入力ミス、または結果の誤った解釈などが発生する可能性があり、実用的な検証が必要である。
一方、最近は扱いやすい表計算ソフトが普及し、煩雑な設計計算から得られた行列データなどを情報通信などを介して容易に表計算ソフト上へ受取ることができるようになったので、数表上での行や列単位の演算がマウス主体で実行できるようになってきている。
本稿は表計算ソフトを利用することにより、高度な固有振動解析に疎遠な技術者でも表計算感覚で機械的に検算ができることに着目した一つの提案である。
検算の基本的な手法には、立体骨組構造の特性を表す変位行列（変位の影響値）を汎用プログラムであらかじめ計算しておき、動的問題を静的問題の延長とすることができる
Stodola の方法 3）と Newmark の方法 4）を組合わせた数値計算法を応用している。
ここでは、SI 単位系で表示された数値計算例を用いた検算表の作成と利用法を主体に考え、計算式が組込まれた検算表などをインターネットを利用して構造技術者に直接提供するものである｡なお、追補として本稿の数値計算法を拡張した固有振動の解析プログラムも計算式、例題とともに提供する。本稿の本文および図は Word97、表および追補（計算式付き）は Excel97 を使用している。ダウンロードされた情報は個人のコンピュタの中で修正追加も含め実際に使用することができる。
関連ファイル「固有振動解析結果の検算（本文）」 Word97「固有振動解析結果の検算（付図）」 Word97「固有振動解析結果の検算（付表）」 Excel97「固有振動解析結果の検算（追補）」 Excel97

１

検算法の考え方

自由振動している任意の骨組構造の運動方程式は周知のように式（１）、式（２）、および式（３）で表される。

1[ M ]{U&&} + [ K ]{U } = {0}（１）固有値、固有ベクトルを介して[ K ]{φ} = ω 2 [ M ]{φ}（２）式（２）が有意の解をもつためには式（３）の固有方程式が得られる。

([ K ] − ω [ M ]) = 02det

ここに

（３）

[M ] ：節点集中質量マトリックス[K ] ：剛性マトリックス{U } ：変位ベクトルω2：固有値（ ω ：固有円振動数）{φ}：固有ベクトル（任意定数に比例）

固有振動解析では式（３）、および式（２）よりｎ組の固有値および固有ベクトルが計算される。
Stodola の方法では式（２）の両辺に1[ F ] をかけ式（４）に書き換える。
ω21{φ } = [F ]⋅ [M ]{φ }ω2
ここに（４）[ F ] ：変位マトリックス（ 変位の影響値、対称）[ F ] = [ K]−1、

[ F ] ⋅ [ K]−1

= [I ]

さらに、式（４）を設計計算等でよく使用される構造物に作用している荷重（自重）と変形の関係に近づけるため両辺にｇ（重力加速度）をかけ、式（５）,式（６）を得る。

g{φ} = [ F ]g[ M ]{φ}ω2（５）

＝ { y}（６）[ ]式（５）右辺の g M は節点に作用している荷重となり（ g をかけているが必ずしも重力[ ][ ] [ ] は節点の変位、さら方向ではない）、また F が変位の影響値であることから F g Mに {φ} をかけた { y} は慣性力に比例した変位とみることができる。
2検算の一つの方法は、結果として得られている数組の ω と {φ} を式(５)に代入し、式(５)が成立つことを確認することである。また、設計計算として実用的に使用される振動2

解析の結果は、一般には１次モード（固有振動数最小、または固有周期最大）を含む数組の固有値と固有ベクトルである。したがって、１次モードが含まれていることが経験的にでも確認されている場合は式(５)の検算のみでよいが、それ以外はｎ組の結果の最初に１次モードが含まれていることを確認しなければならない。
これには式（５）の {φ} を仮定し繰返し計算させるとき必ず１次モードに収束することを利用し数値計算を行うとよい。 Newmark の方法は質量、変位および剛性の分布が２次式[ ]で与えられるものと仮定し、 F を介入させず作表計算により行っている。
本稿では、 Stodola の方法と Newmark の方法を組合わせ、表計算ソフトを利用しやすい表示法で繰返し計算の考え方を述べる。
式（５）、式（６）を書直して、式（７）を、さらに式（７）の左辺と最後の右辺より式（８）、および、式（９）を得る。

g{ ya } = [ F ]g[ M ]{ ya }
ω2= [ F ] ⋅ [ P]{ y a }{ }= [ F ] Py = { y i }yig= const = 2
ω
ya

ω2 =g
( yi / y a )ここに，，n=（７）1−ω2πya '=0ya，

T=（８）

1n（９）

{ ya } ：仮定した変位ベクトル[ P] ：重力マトリックス（対角行列）{ P } ：慣性力ベクトルy{ yi } ： { ya } によって式（７）で計算された変位ベクトルy a ： { y a } の節点ｋの要素
y i ： { y i } の節点ｋの要素
y a ' ： y i を規準化した値

3n ：固有振動数
T ：固有振動周期
ここで、 { y a } として次数に関係なく結果として得られている正しいモードを入力した
場合、繰返し計算しなければ式（７）から式（９）は成立つ。すなわち、次数は確認できないが正しい解の一組であることは確認できる。
一般には { y a } は仮定であるので式（８）は成立しない。そこで { y i } を { y a } と同じ規準値に換算（規準化）し再び式（７）に代入し、式（８）を計算する。これを数回繰返せば式（８）を満足しこのとき１次モードが計算され、式（９）より ω も計算できる。したがって、１次モードのみ、変位マトリックスを利用した固有振動の数値計算をしたことに2なる（高次モードについては４４２

追補で述べる）。

検算の要領

２．１

検算の前提条件

固有振動の解析は、厳密なモデル化を行うものから実用的な近似解に至るまで幅広くある。ここでは解の検算を行うのが目的であるから次の条件で計算を行う。
1)

質量は集中質量とし節点のみに作用する。

2) 変位マトリックス [F]（変位の影響値または単位荷重による変位値）は一般的な骨組構造解析の結果を用いる。したがって、直行座標系において構造特性を示す入力は[F]の計算に用いられる条件に従う。
3) 回転慣性の影響は考えない。したがって[F]は回転変位に関する要素を削除できる。
4) ＳＩ単位系を用いる（重力単位系には記号 kN を tf と読替える）。
２． ２

表計算ソフトを使った検算表の作成と検算手順

検算は式（７）から式（９）を計算することになる。検算表の作成と検算手順を表計算ソフト Excel を使った例として説明する。なお、検算の媒介である[F] は骨組構造の形状（座標）、剛性および境界条件から計算されており、[F]自身に誤りがないことが別の方法で確認されていなければならない。
いま、図―１図―１および表―１
表―１に示すようなバネ支持された単純ばりがあり、表―２表―２のよう図―１表―１表―２な鉛直振動１次、２次および３次モードの解析結果と自重によるたわみ（[F]のチェック用）がある。この例を使用し検算の要領を示す。検算表には数式が組込まれているので以下の説明を演算しながら確認することができる。記号○は説明用に使用する各検算表の見出しである。
１） 任意の Sheet に表―３表―３の変位の影響値[Ｆ]を入力する。
表―３
この例では鉛直変位のみであり、別の解析ソフトによって得られたたわみの影響線（た4わみ図）を Sheet に記入（コピー）する。ここで表―１
表―１のたわみを確認するここによって表―１[Ｆ]が正しいことを検算しておくとよい。この例では[F]*{P}を表計算することになり、表―４（計算式付き）からたわみが一致していることが確認できる。
表―４２）検算表の作成と１次モードの検算
表―５を用意し、式（７）を計算する。
表―５① [Ｆ]と対応する節点記号を記入（コピー）する。
② 仮定したモード ya として表―２表―２の表―２ 1 次モードを入力（コピー）する。
③

表―１の節点集中自重を入力（コピー）する。
表―１④②＊③＝④を計算する（節点１のみ計算しドラッグ）⑤[Ｆ]{Py}={yi}を計算する。この場合行列演算となるので、組込み関数 MMULT（行列の掛け算）を使うと便利である。たとえば、⑤の節点１〜９を選択（ドラッグ）、fx をクリック、MMULT 選択、[Ｆ]の指定（表―表― ３の要素全部をドラッグ、$ 表示にしておくと表全体をコピーしたとき便利）、Py の指定（④全部をドラッグ）、完了

クリック、演算式表示窓をアクティブにして[CTRL]+[SHIFT]+[ENTER]キーを同時に押す。
⑥⑤の各要素（セル）を②で 1.0 とした値に換算（規準化）する。この例では節点５の要素が 1.0 になるようなに⑤の各要素を節点５の要素で割り算する。（絶対番地表示＄を使いドラッグすると便利）⑦

どこかの節点で式（８）(1‑ ⑥／②)を計算し、全体にドラッグする。このとき⑦の要素の絶対値が 0.01 以下（②と⑥が等しい）であれば仮定したモードは実用的に正しいとみなすことができるであろう。（②の要素が 0 または 0 に近い値の節点は考慮しない）

式（９）を計算する。②の要素が 0 または 0 に近い値の節点以外ではどの節点の値を用いてもよいが、振幅の大きそうな節点の方が精度がよい。なお、[F] の表示に係数が考慮6されている場合は式（９）の計算のみに考慮すればよい。この例では 10 である。
ここまでの表計算が表―５表―５であり、モードおよび固有値は図―１
図―１および表―１
表―１に示す構表―５図―１表―１造の解の一組であることが確認できる。
３） 1 次モードであることの確認（確認の必要がないときは省略）
表―５の⑥の数値のみを②に代入する。このときすべての数値が変動しなければ表―５表―５表―５は収束していることになり、1 次モードであることが確認できる。
1 次モードの数値計算を兼ねて繰返し計算を行う場合は表―５表―５の②に仮定したモードを表―５代入し繰返し計算を行う。説明のため表―５
表―５を表―６表―５ 表―６にコピーして（$表示注意）繰返し
表―６計算の要領をのべる。
表―６の②に振幅の大きそうな節点の値を1.0 とした適当な数値を代入する。ここで表―６は節点５を 1.0 、他はすべて 0（空白）とした。
仮定値②が正しければ式（８）により⑦はすべて 0 になる。この例ではまだばらつきがある（⑥と②が等しくない）ので、⑥全部の数値のみ（数式を除く）を②に代入する。

5たとえば、⑥全体をドラッグ、コピー、②の先頭を指定（クリック）、形式を選択して貼り付け、値、OK クリック。この操作によって新しい仮定値②による⑦の値が同じ表の中で計算される。形式を選択して貼り付け以降を数回（この例の場合３回）繰返せば（繰返しボタン使用）⑦はすべて 0 となり１次モードに収束する。繰返し数を増せばいくらでも精度は向上するが１次モードで収束させるには⑦は 0.001 以下で十分であろう。このとき、固有値も計算され、表―３の表―３の解析結果と一致していることが確認できる。
表―３の
４）２次モードの検算
表―５または表―６
表―６を数式ともコピーし表―７表―７（規準節点が変わるため⑥と固有値の表―５表―６表―７数式を修正）を用意する。
② 仮定したモード ya として表―２表―２の２次モードを入力（コピー）する。
表―２⑦は 0.01 以下（②と⑥が等しい）となる。
固有値も計算され、表―３の表―３の解析結果と一致していることが確認できる。
表―３の
５）３次モードの検算
表―５または表―６表―５表―６を数式ともコピーし表―８表―８（規準節点が変わるため⑥と固有値の表―６表―８数式を修正）を用意する。
② 仮定したモード ya として表―２表―２の３次モードを入力（コピー）する。
表―２⑦は 0.01 以下（②と⑥が等しい）となる。
固有値も計算され、表―３の表―３の解析結果と一致していることが確認できる。
表―３の
２． ３

検算表の作成と検算手順のまとめ

２．２では簡単な構造で説明したが、一般の平面構造または立体構造に対しても若干の２．２修正を加えればよい。たとえば、図−１
図−１は振動モードの成分が鉛直方向のみであるが、一図−１般には２（３）次元である。よって次のように拡張すればよい。
１）構造形式に関係なく検算表のフォームは一定であるので表―５
表―５（本稿で使用されてい表―５る例題のどの検算表でもよい）を数式ともコピーする
２）検算対象の変位節点数に合わせて、検算表の行を増減する（列は変わらない）。すな
わち検算表が縦長になるのみである。ここで、変位節点数とは節点の変位成分数の総和である。
３）①に検算対象の節点記号を入力（コピー）する。[F]の行または列の並び順と一致していなければならない。また、[F]の行および列の要素がすべて 0 の節点（支点）は検算表の行、[F]の行および列を削除してもよい（モード値は 0 ）。
４）②に基準値を 1.0 とした検算対象のモード値を入力（コピー）する。基準値はすべての成分の中で変位の大きそうな節点の値をとるのがよい。
５）③は変位成分に関係なく１つの節点について同じ自重（本来質量であるがここでは自重としてる）を入力する。なお、[F] で考慮された成分であっても振動モードとして無視される成分がある場合（慣性力④が作用しない）は節点記号の③は 0 とする。
６）[F]の領域の変更に伴う⑤の数式の修正する。２．２
２．２6１）⑤参照。７）規準節点の変更に伴う⑥の数式の修正する。２．２
２．２１）⑥参照。

８）検算表の表示桁数、枠線などは検算対象ごとに調整する。

３検算例

３．１

連続補剛桁を有するアーチ橋
連続補剛桁を有するアーチ橋

図―２は連続補剛桁を有するアーチ橋（仮想）である。表―９表―９の解析条件に対して表―表―
図―２表―９10 の解析結果が得られている。
表―11
図―２、表―９を計算し Sheet に記入したものである（水平変位成
表―11 は図―２
図―２ 表―９から[F]表―９分 U、垂直変位成分 V）。表―
表―12表
表―12 は[F]自身の検算として鉛直変位 [F]{P’}を表計算し、表―９のチェック用数値と照合して座標、剛性および支点条件等に誤りのないことが確認さ―９れている｡
表―13
表―15
表―９の１次、２次およ
表―13、表―13 表―14表―14 および表―
表―15 は検算表であり､それぞれ②に表―９表―９び３次のモ−ドを入力（コピ−）し、⑦が 0 になることを確認している｡
とくに表―
表―13表―13 では⑥の数値を②にコピ−しても数値が変動しないので１次モ−ドに収束していることが確認できる｡なお、表―
表―13表―13 を用いて１次モ−ドの数値計算をする場合には、たとえば②の V13=1.0、V21=‑1.0、他は 0 （空白）を初期値として繰返し計算を行えば 6 回程度で収束し表―
表―13表―13 と一致する。

３．２

曲線格子桁

図―３は多角形で近似したは曲線格子桁（仮想）である。表―
表―16図―３表―16 の解析条件に対して
表―17表―17 の解析結果が得られている。
表―18図―３、表―
表―19表―18 は図―３図―３ 表―16表―16 から[F] を計算し Sheet に記入したものである。表―
表―19 は[F]自身の検算として鉛直変位 [F]{P’}を表計算し、表―
表―17表―17 のチェック用数値と照合して座標、剛性および支点条件等に誤りのないことが確認されている｡
表―20表―22表―17表―20、表―
20 表―21
表―21 および表―
表―22 は検算表であり､それぞれ②に表―
表―17 の１次、２次および３次のモ−ドを入力（コピ−）し、⑦が 0 になることを確認している｡
とくに表―
表―20表―20 では⑥の数値を②にコピ−しても数値が変動しないので１次モ−ドに収束していることが確認できる｡なお、表―
表―20
表―20 を用いて１次モ−ドの数値計算をする場合には、たとえば②のδz-14 = 1.0 、 他は 0（空白）を初期値として繰返し計算を行えば３回程度で収束し表―
表―20表―20 と一致する。

３．３立体アーチ

図―４は簡単な立体アーチである。表―表―23表―24図―４表―23 の解析条件に対して表―表―24 の解析結果が得られている。
表―25図―４、表―
表―25 は図―４図―４ 表―23表―23 から[F] を計算し Sheet に記入（節点１と７はすべての変位7が 0 であるので削除）したものである。表―
表―27
表―27 は[F]自身の検算として自重による鉛直変位 [F]{P’}を表計算し、表―
表―23表―23 のチェック用数値と照合して座標、剛性および支点条件等に誤りのないことが確認されている｡
表―28表―30表―28、表―28 表―29表―29 および表―
表―30 は節点１と７を削除（0 として考慮してもよい）した検算表であり､それぞれ②に表―
表―24表―24 の１次、２次および３次のモ−ドを入力（コピ−）し、⑦が 0 になることを確認している｡
とくに表―表―28表―28 では⑥の数値を②にコピ−しても数値が変動しないので１次モ−ドに収束していることが確認できる｡なお、表―
表―28表―28 を用いて１次モ−ドの数値計算をする場合には、たとえば②のｙ4 = 1.0 、 他は 0（空白）を初期値として繰返し計算を行えば８回程度で収束し表―
表―28表―28 と一致する。
変位の種類が多くなる場合では解析結果として表―
表―23
表―23 と節点の順序が異なった表示の場合がある。この場合は、[F] の表示を表―
表―26表―28表―26 のように並べ替え、表―表―28、表―28 表―29表―29 および表―
表―30表―30 と同じように計算することができる。検算表と[F]の節点記号の並びが同じであればどのような並びでも検算できる。表―
表―31表―26表―31 は表―
表―26 の[F]に対する検算表である。

４追補４．１４．１

固有振動の計算

１ において式（７）および式（８）は 1 次モードのみに収束する。これを 1 次モード以外でも収束するように拡張するには、掃出しマトリックスを用いて変位マトリックス[F] を更新する方法がある。いま、n+1 次で収束する[F]を計算する場合次の式となる。

[ ff n ] = [ ff n −1 ] −1{P}{φ n }{φ n }T
Mp n（11）

Mp n = ¦ φ ni ⋅ P（12）[ Fn +1 ] = [ F ] ⋅ [ ff n ]

（13）

2

ここに

[ ff n ] ： ｎ次までの変位成分を除去した掃出しマトリックス
ただし、 [ ff 0 ] = [ I ]（単位行列）

Mp n : 規準化係数、ｎ次モードの２乗とＰの積の変位節点の総和[ Fn +1 ] ： n+1 次で収束する変位マトリックス
式（13）が計算できれば、式（７）から式（９）の[F] と読替えることによって収束させることができる。すなわち、1 次から順に式（７）から式（13）を繰返し計算し任意の次数まで計算することができる。

8４．２
４．２ 表計算ソフトによる固有振動の計算
式（７）から式（13）を Excel の表計算で直接対応するには、データの扱い量が多いのでやや問題がある。ここでは簡単なプログラムの知識が必要であるが、参考として Excelのマクロで記述した検算表を提供する。
マクロのリスト、プログラム仕様および計算例は関連ファイル「固有振動解析結果の検算（追補）」を参照されたい。
あとがき
設計計算と直接関係のない表計算ソフトを使い、一例として Excel による固有振動解析結果の検算法を述べた。
表計算ソフトを有効に使うには、すべてのデータが Sheet に記録されていなければならない。データ数が少なければ直接手入力することもできるが、一般には他のプログラムで計算された[F] や固有モードなどを手入力することは困難である。本稿では Excel のデータに変換して表計算に取込むことを前提にしている。
なお、追補としてまとめた固有振動の数値計算は、コメント行を多くとった Excel マクロ (VBA)で書かれたファイルを直接提供するので、ダウンロードしたファイルの修正、追加も含め利用していただきたい。
最後に、本稿が ISO9001 の要求項目の一つである

設計検証の一部にでもお役に立てば幸いに思うところである。

参考文献
1) チェモシェンコ（谷下、渡辺訳）：工業振動学、東京図書、1956.5
2) 小坪清真：土木振動学、森北出版、1973.4
3) Clough,R.W.andPenzien,J. （大崎、渡部、片山訳）：構造物の動的解析、科学技術出版社、1978.5
4) 成岡昌夫：ニューマークの数値計算法 梁・柱の曲げ、振動、座屈に関連して、技報堂、1978.65) 中平、夏秋、小澤：Newmark の数値計算法によるはりの固有振動周期の計算、土木会
誌、1983.6,pp.48‑54
6) 中平、小澤、成岡：四則演算のみによるトラスの固有振動数の計算、土木学会誌、
1986.3,pp.29‑367) 西岡隆：土木工学基礎シリーズ４ 構造振動解析、培風館、1987.48) 中平、成岡、小津：ニューマークの数値計算法による梁・トラスの固有振動数の計算、橋梁と基礎、1988.2,pp.11‑169) 安藤、織田：骨組構造の固有振動解析結果の検算、技報たきがみ、2001.5,pp71‑809