フジタ技術研究報告第 50 号 2014 年

有限要素法による床振動解析システム「揺れイザーⅢ」の開発

仲沢近藤

概武志貴士

増田圭司

*1

要

構造物の大型化や建物の内部空間を有効に利用することなどの要請から長スパンの梁や部材の軽量化が必要とされている。
その際、構造的な強度だけでなく常時発生する振動が問題となる場合もある。これに対応すると同時に、コストダウンを目的として部材に鉄骨を用い、端部を鉄筋コンクリートで支えるハイブリッド構造梁が使われるようになってきている。
一方、常時振動問題の予測手法として「揺れイザー」を開発しているが、ハイブリッド構造は複雑な部材であるので、「揺れイザー」での計算手法であるレイリー・リッツ法では適用に限界が生じる。また、床スラブの開口に対応できないことや、床厚さが均一でない場合にも対応できないなどの問題も生じていた。そこで、均一部材だけでなく部材内で特性が変化する場合にも適用できるような重みつき残差法による解析手法を開発した。ここでは、その内容を固有値解析を対象として記述する。
なお、本システムは数値解析分野で使われるメッシュジェネレータを使用せず、実務で馴染みの有る２次元図面のイメージで３次元データを作成できるプリおよびポストプロセッサーを装備している。
本文はその概要と適用例を示すものである。

Development of finite element analysis system for floor vibration named “Yurayzer Ⅲ”

Abstract. Recently, there has been increasing demand for buildings with maximum usable space. As a result,
there is a tendency towards construction using long spans and lightweight members. This leads tochallenges both in structural strength and environmental performance, including issues such as floor
vibration. To address this problem and simultaneously reduce construction costs, hybrid members
combining steel beams with reinforced concrete end sections have been used for constructing
buildings.
The authors have previously developed a prediction analysis method of floor vibration problems
named “YureizaⅠ&Ⅱ” that is based on Rayleigh-Ritz's method. However, the calculation method islimited to members with uniform mechanical characteristics. So we developed another calculation
method in this study. That method is called Galerkin method based on weighted residual procedure.
This system has also functions to provide input data for calculations and to show some graphics ofcalculation results. Especially, users can easily provide 3 dimensional input data for calculations
based on images of usual 2 dimensional drawing by using this system.
In this paper describe the brief summary and example of application.

キーワード： 床振動、有限要素法、固有振動数*1 大和ハウス工業株式会社 総合技術研究所

振動モード－1－

フジタ技術研究報告 第 50 号

§１．はじめに構造物の大型化や建物の内部空間を有効に利用することなどの要請から長スパンの梁や部材の軽量化が必要とされている。また、コストダウンを目的として部材に鉄骨を用い、端部を鉄筋コンクリートで支えるハイブリッド構造梁も使われるようになってきている。
構造物の建設には、必要強度を満たすことは当然として、居住性などの検討も必要となる。
フジタは、この居住性の検討に対して従来から「揺れイザー」と称する床振動解析システムを開発してきた。この解析
では、レイリー・リッツ法に基づいて、基本的には解析対象

 12 EI xx
 l30
00 6 EIxx
2 l0
 12 EIxx
l3
0

0

0
 6 EI xx
 l2
0


00

12 EI yy0l3l2EAl0
6 EI yy00
4 EI yy

6 EI xx
l200000000

4 EI xx
l00000000

12 EI yy0l3
0
6 EI yy

EAl0



6 EI yyl202 EI yy00l300

6 EI yy0l2EA
l0l26 EI xx
l206 EI yy002 EI yy0l000GK
l000006 EI
 2 xxl012 EI xx
l306 EI 2 xxl00000000002 EI xx
l00000l012 EI yy

6 EI xx
l20l2

12 EI xx
l300l



00l2

0
6 EI yyGKl00l3000
12 EI yyEA
l6 EI yyl202 EI xx
l6 EI yyl204 EI yy
l0006 EI xx
l20004 EI xx
l00000
0 
0 
0 
0 
0 
GK 
l 
0 
0 
0 
0 
0 GK l 

図 1 線材の剛性マトリックスが均一な部材に限定される。つまり、上記のようなハイブリッド構造にそのままの形で適用することは困難となる。

z

そこで、ハイブリッド構造にも適用できるようにガラーキンy法による有限要素法での床振動解析プログラムを開発した。
さらに、独自のプリポストも開発し、両者を連携させることで利便性の向上も考えたシステムとなっている。
本文では、その内容を例題とともに記述する。

§２．計算手法
計算には、面材と線材を用いる。梁や柱は線材としてモxデル化し、ビーム要素で表現する。床は面材でモデル化し、薄板理論に基づく三角形要素を用いる。なお、この解析で

図 2 三角形要素の基準座標は、通常の変位３成分と曲げ３成分から６自由度を各節点は持つこととなる。のような面外変形による曲げを生じる。ここでは、面内変形

以下に定式化を示す。

に対して定ひずみ三角形要素を用い、面外変形においては面積座標を用いた計算方法を採用する。

2.1 線材 1)

三角形は３節点で構成されるが、線材と同様に１節点６自線材としてのビーム要素には、従来から提案されている由度であるので、要素としては１８自由度を持つこととなる。

古典的なものを用いる。すなわち、軸力と同時に曲げをビ

ここでは、図2のように、 xy 面内に三角形要素が横たわり、ーム要素は負担する。その剛性マトリクスを図１に示す。ここたわみは z 軸方向を基準に生じる状態を基準とする。に、 E は弾性係数、 I は断面２次モーメント、 G はせん断

各節点の自由度には、 (xy ) 面内に生じる面内変位(u , v) と (z ) 軸方向に生じるたわみ (w) および各軸を基準としたたわみ角 ( x ,  y ,  z ) がある。以下に、これらを分弾性係数、 K はサン・ブナンのねじりこわさ、 l は要素の長さをそれぞれ示している。
ビーム要素は２節点で構成されることから要素としては１けた形で計算に用いる離散化を記述する。なお、薄板理論２自由度を持つこととなる。また、図１は Z 軸方向を材軸にでは  z が明示的に考慮されていないことに注意を要する。

した表現であるが、材軸がこれ以外の向きの場合、対応し① 面内変形 2)た座標変換を施す必要が生じる。なお、ここでは剛性のみ

面内変形には、定ひずみ三角形要素を用いる。この場合、記しているが、質量マトリクスについては文献 1)を参照され要素内の変位分布を表現するために導入する形状関数はたい。

以下で与えられている。

u  1   2 x   3 y(1)v  1   2 x   3 y式(1)の未知係数  i (i  1～3）を変位成分 (u , v ) と座標( x, y ) で表し、未知変数を消去する。このとき用いる座標

2.2 面材
面材には構造全体形状の自由度に対応できるために三角形要素を採用する。面材には、面内変形とともにたわみ－2－有限要素法による床振動解析システム「揺れイザーⅢ」の開発

は三角形要素を構成する各節点座標となる。変位成分 (u )について考えると、(2)ui   1   2 xi   3 yi (i  1～3）となり、変数 3 個に対して式が 3 本立ち  i (i  1～3）を確定的に得ることが可能となる。この手順は変位成分 (v) についても同様に適用される。計算された  i (i  1～3）を再度式(1)に代入することで、未知変数が消去された形状関数を得ることができる。

図３ 三角形要素の面積座標得られた形状関数を次式のように書く。

{ }  B{u}{u i }  N m ( x, y ){u mp }(3){vi }  N m ( x, y ){vmp }pp
ここに、 (u m , vm ) は離散化された後の各節点の変位を表式(5)を組み込んで式(4)を書き換えると、通常得られている要素剛性マトリクスが以下のように計算される。

す。 以上で、面内変形に対する形状関数を確定できたこ

 { f }ds   B DBdv{u }Tpとになる。

spv(6)

式(6)を全要素重ね合わせたものが面内変形における連立

次に、変位を求めるために場の方程式について考える。
これに使われる概念は、次式に示す「仮想仕事の原理」や１次方程式となる。

「クラペイロンの定理」といわれるものであり、つり合い条件② 面外変形 3)を仕事やエネルギーで表現したものである。

面外変形における目的変数は、薄板理論の仮定から

TTs( z ) を除いた ( w,  x ,  y ) となる。各節点で面外変形は(4)

 {u} { f }ds   { } { }dv

3 自由度となることから、三角形要素の場合各要素 9 自由度v{ f }  { f x , f y } ：外力ベクトル、 {u}  {u, v} ：変位ベクトル、 { } ：応力ベクトル、 { } ：ひずみベクトル、 V ：体積、 S ：表面積、  ：増分あるいは変分TTもつこととなる。
面積座標は図 3 に示す座標系であるが、それを使った形状関数は以下のように導入される。

ただし、外力ベクトルと変位ベクトルは全要素構成節点分w  {}T { }考慮する。

ここに、 { } ：面内変形と同様な未知係数

応力ひずみ関係は、次式のように得られている。

{} は面積座標を使って次式のように導入される。

{ }  D { }

(5)

1
12 1 2{ }    2 3  

3 1
 1 22   2 12 22
 2 3   3 2 
 3 12   1 32 { }T  { x ,  y , xy } { }T  { x ,  y ,  xy }ここに、

なお、今の場合は平面応力が仮定されるので、応力ひずみマトリクス

D  は以下を用いることとなる。

10 
D  E 2  1 0 1  1  0 0
2 
ここに、 E ：弾性係数、 ：ポアソン比
ここに、式(4)のひずみベクトルは形状関数を各成分の定義どおりに微分することで得られる。

mA

m1n2

p
3(8)

m!n! p!

    dA  2 A (m  n  p  2)!

N
Nu mpy v mpxy
N m p  N m p
 xy um vmx
yx m(7)(9)

以下の手続きで未知係数{ } を求める。

最終的に求めたいのは ( x, y ) 座標における剛性であるが、一旦式(7)の未知変数を面積座標で求めて ( x, y ) 座標

 

これらは、いわゆる B マトリクスを導入することで、以下の系に変換する手順で計算する。節点の自由度はたわみとたようにまとめて表現される。

わみ角であるので、まず各節点のたわみ角を面積座標によ－3－フジタ技術研究報告 第 50 号

   Bd  とおく。面材のひずみベクトルは、式(15)を って微分し、自由度にはめこむ。

 T  w, xx wi (d  ) i   wi , 1 w ,  i 2

i  1～3(10)であるので、面内変形と同様に形状関数を微分することでこれらは次のように計算される。

   P Bd 

ここにカンマは対応する変数での微分を表す。
面積座標で得られた自由度と ( x, y ) 座標系での自由度のここに、関係は、ヤコビアンマトリクスによって次式のように得られて     x   1   J     y   2 いる。

 , yy2 , xy   D   DP Bd   xJ    ( x, y )   x1
 ( 1 ,  2 )   2P T   , xx(17)よって応力ベクトルは次のように得られる。

(11)

ここに、 J はヤコビアンマトリクスで次式のように計算される。

2 w, xy w, yy

y  1 y 
 2 

(18)

ここに、応力ベクトルは曲げモーメントで T  M xMyMZ

材料定数マトリクスは
(12)

0 
1 D  D  1 0 
1  0 0
2 
3EhD12(1  2 )これより、

1 0 (d  ) i   (d ) i
0 J 
(13)
式(13)のような形で面積座標と ( x, y ) 座標での自由度の関係が得られる。ただし、式(13)は３×３の係数行列を簡易表記していることに注意を要す。

ここに、 h は板厚でその他は面内変形と同じものを示して一方、形状関数式(7)を節点自由度式(10)に代入したものいる。

を、全要素構成節点について、式(17)および式(18)を式(16)に代入することによって、面材 ( A )1 ( d )  ( A ) 2    A  
( A ) 3 に対する要素剛性マトリクスは次式のように得ることができる。
(14) f   A BT P  DP BdAd T(19)

なお、式（１９）の実際の積分は Gauss の９点積分で実施した。さらに詳しい記述は文献 3)に譲る。と表す。ここに、 A i  i , 1 i ,T 2 Ti
T

以上、面材においては面内変形 (u , v) と面外での変形( w,  x ,  y ) および薄板理論では明示的に考慮しない( z ) が適切な自由度で重ね合わせされ各要素で 18 自由度をもつ剛性マトリクスが作成される。



これより、未知係数  は次式のように求めることができる。

さらに、線材の剛性と面材の剛性を適切に重ね合わせて構造物全体の剛性マトリクスを作成し、境界条件の処理を

1 0   A  d   A  d  (15)
0 J 11

施して目的とする計算を実施し、解を得ることができる。
なお面材における質量マトリックスは単位体積重量

このようにして式(7)の未知係数が計算される。

を  とし、式(1)や式(7)の形状関数を使って次式で計次に要素剛性マトリクスを求める。面材に対する場の方程算される。

式は次式のように得られている。

M     N N dVmnmnVd T  f    A  T  dA(16)上式は必要に応じて集中質量マトリクスに縮約することも可能であるので、本システムではこれを採用している。

－4－

有限要素法による床振動解析システム「揺れイザーⅢ」の開発(a) 側面図(b)平面図図４ 解析対象

図 5(a) 1 次モード(case1)

図 6(a) 1 次モード(case2)

図 5(b) 2 次モード(case1)

図 6(b) 2 次モード(case2)§３．適用例

図 5(c) 3 次モード(case1)

図 6(c) 3 次モード(case2)

計算対象の概略を図４に示す。図は対象建物の３階部分であり、そのうち 1 スパンのスラブのみ対象としたもの（case1）と周辺のスラブを考慮したもの(case2)とを計算した。

前述の計算手法の適用として固有値問題を考える。対象構造物が大規模で挙動を簡易に知りたい場合、その一部を切り出して計算することは比較的自然な発想であるが、切り3.2 計算結果出した部分の境界の影響を調べる必要がある。ここでは、そ図 5 は case1 に対して得られた固有モードの 1 次から３の境界の影響を調べることとする。なお、ここでの固有値解次までの計算結果を表す。また、図６は case2 に対して同じ析の計算にはハウスホルダー法を使用している。

く得られた１次から３次の固有モードを表す。
図 6 に示す固有モードに関しては、3 次モードの形状に相違が生じているように思われるが、振動特性として卓越す

3.1 固有モードに関する計算－5－フジタ技術研究報告 第 50 号

いて実測された値と計算との対比を図 8 に示す。さらに、これらの統計的比較を表１に示す。両者ともに変動係数は10％程度に納まっており、良好な精度で計算されているものと思われる。

§４．おわりに
本文においては開発した「揺れイザーⅢ」内部での計算概要ならびに得られた結果を示した。「揺れイザーⅢ」の特徴はデータの作成や結果の表示にいわゆる「メッシュジェ

図 7 固有振動数に対する解析対象ネレータ」を使うことなく簡易なインターフェースで実施できることに特徴がある。また、若干の例題を示しているが、その結果も良好であるので、今後は実物件での設計や検討に適用したいと考えている。
参考文献
1) 川井忠彦、藤谷義信：振動および応答解析入門、有限要素法の基礎と応用シリーズ５、培風館
2) 例えば、O.C ツェンキーウィッツ著、吉識雅夫、山田嘉昭監
訳：マトリックス有限要素法三訂版、培風館
3) 山田嘉昭編：マトリックス法の応用、東京大学出版会

図 8 固有振動数の比較（実測と計算）

ひ と こ とフジタでは以前から床振動

表 1 固有振動数の統計的結果の問題に多くの先輩達が携わっている。今回、あらためて解析システムを構築して、ハイブリッド構造のように複雑な構造部材に適用範囲が仲沢 武志広がった。今後活用される事を願う。

る 1 次や 2 次モードに関してはある程度同様な結果が得られている。解析の目的や重要度にもよるが、1 スパンの解析でも簡易的な検討としては有効な場合もあると考えられる。
3.3 固有振動数に関する計算
固有振動数に関しても図７に示すように解析領域を変えた計算を実施している。すなわち、中央１×1 スパンのみの計算と周辺のスパンを考慮した 3×3 スパンの 2 ケースを設定した。
計測自体の内容は別の報告に譲るが、固有振動数にお－6－