第 6 章 有限要素法 ( その 2) 振動問題を有限要素法で解いてみよう. 振動方程式は式 (3.35) で与えられ (6.1) [ K] { d ( )} + [ M] d ( t ) { } F( t ) t = { } そのときの質量行列は式 (3.32) で T M N N d V (6.

Transcription

1 第 6 章 有限要素法 ( その ) 振動問題を有限要素法で解いてみよう. 振動方程式は式 (.5) で与えられ (6.) [ K] { d ( )} + [ M] d ( t ) { } F( t ) t { } そのときの質量行列は式 (.) で M N N d V (6.) である. 6. 固有値解析法 [ ] ò { } { } r 系が調和振動をしている場合, 外力はなく式 (6.) は K d t + M d t (6.) でそのとき節点変位 d は と表せるので加速度は である. 従って式 (6.) は すなわち となる. これを書き直すと { } [ ] { ( )} [ ] ( ) wt { ( t )} { d } e d (6.) { wt ( t )} w { d } e d [ K] { d } w [ M ] { d } { } [ ] { d } w [ ] { } K M (6.5) d [ K ] [ M]{ d } [ I]{ d } w (6.6) / w としている理由は, 振動問題では実際には となって標準的な固有値問題となる. 固有値をわざわざ低い振動数が問題となるのに対して多くの固有値解析のプログラムが固有値の大きい方から求めるように作られているため, 逆数をとってそのプログラムに対応させるためである. 6.. 固有振動モードの直交性 式 (6.6) を標準形 [ ]{ } l[ I ]{ } と書いて, 固有ベクトル の直交性を考察する. 第 番目の固有値を A (6.7) l, 固有ベクトルを { } として [ A ]{ } l { } 第 番目については 6

2 これら つの式の上の式に左から { [ A ] { } l { } }, 下の式に { } をかけて引き算すると ( ) { } [ A]{ } { } [ A]{ } ( l l ) { } { } [A] が対称であれば ( 対称行列の固有値と固有ベクトルは実数 ) 上の式は ( l l )( { } { } ) となり, l ¹ l であるので ¹ のとき { } { } となる. これを固有ベクトルの直交性という. さて, 一般の振動方程式 (6.5) について考えると [ ]{ } l [ M]{ } (6.8) K (6.9) 第 固有値について [ K ]{ } l [ M ]{ } 先ほどと同様に左からそれぞれ { となる. すなわち },{ } ( l ) { } [ M] { } l ¹ のとき { } [ M] { } となる. のとき固有ベクトル { } あるいは ( ) をかけて引き算をして (6.a) { } [ M] { } [ I ] の絶対値を操作して (6.b) { } [ M] { } [ ] (6.c) と対角行列化することができる. これを質量行列に関する直交性という. 式 (6.9 をかければ { } [ K] { } l { } M] { } ) の前から { } (6.) となるので [K] に関しても直交性が成立する. 6.. 境界条件がない場合 ( 自由 自由 ) の振動解析 拘束がない場合には [K] の逆行列がとれない. これは物理的に当然のことで剛体運動が含まれているわけである. 静的な構造物の場合, 荷重がまったくとれないので解析以前の問題で拘束条件を課せるような構造に改造すべきであるが, 宇宙空間に浮かんでいる衛星のように拘束を受けない構造の振動は特別に考える必要がある. 自由 自由な振動解析の技法として次のような方法がある. 整合質量行列は正定であるので逆行列が計算できて w M d w M (6.) なる恒等式を式 (6.5) に加えて [ ] { } [ ] { } d ([ K] w [ ] ) { d } ( w w ) [ M]{ } M (6.) d として 6

3 ([ K] w [ M] ) { d } I ( ) [ ] { d } w w (6.) なる固有値問題を解くことにする. 固有値は w に近い値から計算されることになる. 式 (6.) を加えたことにより固有ベクトルは影響されないことは式 (6.) をみれば明らかであろう. この技法をシフティング (shftng) と呼ぶこともある. 6. 次元骨組構造の質量行列 6.. 梁の曲げの質量行列 梁の曲げの変位関数 (.5) { }{ d} w N (.5) を質量行列の定義式 (.) に代入して [ M ] ò { N( )} { N ( )} rad (6.5) を得る. ra 質量行列に関しては直感的に質量を節点に等分に振り分けて [ ] (6.6) ra M (6.7) とする集中質量行列 (luped ass atr) も用いられることが多い. 本来, 第, 対角要素は であるが, 質量行列を正定とするため適当な小さな値 ( この場合,(6.6) 式の (,) 成分と (,) 成分を加えた ) を入れることもある. これは前節のような自由 自由の境界条件の解析で集中質量行列を使う場合に, この集中質量行列を正定行列にしてシフティングの技法を使えるようにするためである. この行列を採用することの利点は対角行列となるので計算機における記憶容量と演算回数が少なくてすむことである. 式 (.) で与えられる質量行列は剛性行列と同じ変位関数から得られるので整合質量行列 (consstent ass atr) と呼ばれる. 両者の計算精度が気になるところであるが, 要素分割を十分におこなえば両方とも正解に収束することが保証されている. 特に整合行列を使用した場合, 振動数が高い方から正解に一様収束することが知られている. 6.. 例題 : 片持ち梁の振動 図のような一様梁の曲げ振動を考える. 変位関数として静的な関数を用いているため, 梁に分布している質量を節点上での値に近似するためある程度の要素分割をしなければならない. 今は先の例題と同じく 要素に分割する. 6

4 この解析を付録のプログラム BEAMD により解析して式 (5.) と比較すれば下表のようになる. ただし梁は 6 ページにあるような長さ, 直径 c のアルミの円断面梁で 9, A., E 7 Kg / s, r 8Kg /, 8 I. 785 である. 式 (5.) での l の値は低い方から.857,.69, である. 高次の振動モードが必要な場合にはもっと分割しなければならない. 次数 BEAMD (5.) 式 図 6. 片持ち梁の固有振動数の比較 6.. 縦振動と捩り振動の質量行列 である. 棒の縦振動の質量行列は質量を / づつふりわけて / [ ] ra ねじりの質量行列は極 次モーメント p である. 5 M (6.8) / / [ ] ri 図 6.: 片持ち梁 I に密度 r をかけたものを当分して M P (6.9) / 6.. 次元での質量行列 次元では剛性行列を作ったときと同じようにして式 (6.6),(6.8),(6.9) をまとめて 6

5 65 [ ] A I r / A M p A I / p (6.) となる. これに対応する節点変位は (.7) で表される. 全体剛性行列に組み立てるときには式 (.) と同様に [ ] [ ] [ ] [ ] M M (6.) とする. 6. 大規模固有値計算法通常の計算ライブラリーに付属している固有値計算法はヤコビ法やハウスホルダー法であり, これではせいぜい数 自由度の固有値しか計算できない. 有限要素法の場合, 解析自由度は数百から数万であるがすべての固有値を計算する必要はなく, そのうちの低い値から数十程度でよいので, その特性を活かした固有値計算法が有限要素法の発展とともに開発されてきた. それらは べき乗法 ( パワー法 ) サブスペース法 二分法であり, それぞれ特徴があるが, 最も強力でよく使われている方法がサブスペース法 (subspace ethod,

6 または同時反復法ともいう ) である. 6.. べき乗法 べき乗法 ( ベキはパワーなのでパワー法ともいわれる ) は固有ベクトルを求めるときの反復解法で, 固有値の大きいほうからつの固有ベクトルが得られる. このため座屈解析における固有値解析に採用されることが多いがもちろん振動問題にもよく使われる. 他の固有値解法に比べて丈夫さ ( 条件を選ばない. たとえば対称行列でなくてもよい, 複素問題でもよい, 重複固有値があってもよいなど ) がとりえである. 振動方程式 (6.6) を (6.7) のように表して [ ] { } l{ } A (6.7) とする. ここで l の大きい方から求められれば振動数 w の小さい方から計算できることになる. ということで, 振動問題に関していえば逆べき乗法と呼ぶことができる. u とすると, 任意のベクトル { } は さて, 算法であるが, 正しい固有ベクトルを { } この両辺に [A] を左からかけると となる. この操作を 回行うと A C l 左辺は { } C { u } + C { u } + + { } C n u n [ A ]{ } C [ A]{ u } + C [ A] { u } + + C n [ A] { u n } C l { u } + C l { u } + + C l { } [ ] { } { u } + C l { u } + C l { u } [ A ] { } l { } であるので æ l ö æ l ö ç ç è l ø è l ø を得る. l を一番大きな l と同じ値にすると æ l ö ç << for ³ è l ø であるので { } C { u } n { } C { u } + C { u } + + C { u } として { } は 番目の固有ベクトル { } u に収束する. EA [ K ] [ M] n n u n n n æ l ç è l ra n ø ö n n 5 66

7 EA r A EA r A w M K [ ] EA EA K [ ] [ ] 5. E M K A] [ r [ ] í í í í í A E r w E r w r w E.9 r p p w E. f 9

8 p. 56 E l l r 厳密解 固有ベクトルが得られれば固有値は式 (6.7) の左から { } { } { A}{ } l { } { } をかけることにより として得られる. この方法では第 固有ベクトルしか求められないが, 第, 第 も抜き取り法 ( デフレーション ) により求めることができる. これは固有ベクトルの直交性を利用するもので固有ベクトルを u しておいて に正規化 { } ( ) A A l u u 対称行列のとき 新 [ ] [ ] { }{ } A A l u v 非対称行列のとき 新 [ ] [ ] { }{ } として 番目以降の固有ベクトル計算を行う. ここに { v} は [ ] A に対する固有ベクトルである. ただし, この抜き取りの操作を重ねていくにしたがって誤差が積み重なるので多くの固有値を計算したいときは次のサブスペース法の方がよい. 6.. サブスペース法 サブスペース法は初期固有ベクトルを 個用意する. はせいぜい数 で行列の次元 n( 数 数 程度 ) に比べて著しく小さい値である. サブスペース法では n n の固有値問題を の小さな行列に変換して解く方法である. まずこの初期固有ベクトル { }(,,..) を並べて行列 [X] を作る. すなわち [ X] [{ }{ }{ } { }] である. 次に [Z][M][X] ををつくって [K][Y][Z] を解いて [Y] を求める. すなわち [ Y] [ K] [ M][ X] を得ることになるが, これは先ほどのべき乗法における [A] を左から掛けていく操作と同じことであることに注意されたい. この Y を使って M と K を の小さな行列に変換する. すなわち としてこの固有値問題 [ K ~ ] [ Y] [ K][ Y],[ M] [ Y] [ M][ Y] [ K ~ ] { y} l[ M ~ ] { y} をヤコビ法などで解く. この固有ベクトルをまとめて [ P] [{ y }{ y }{ y } { }] y をつくりもとの空間へ変換する. すなわち [X][Y][P] としてこの手順を繰り返す. ベキ乗法と同じ原理により収束していくので固有ベクトルで収束性を判定してもよいし, ヤコビ法の計算で出てきた固有値で判定してもよい. 詳しくはこの方法の開発者の Bathe の教科書を参照されたい. 68

9 6.. 縮約 縮約とは行列の大きさ ( 次元 ) を小さくするテクニックで要素行列の段階での, また, 解析前の全体行列の段階でもしばしば使われる方法で, それらの中で Guyan の縮約法 (Guyan's Reducton) が簡単でもっとも広く使われている. この方法は静的な剛性方程式 (.) において [ K ]{ d } { F} において残しておきたい節点変位 (aster) をd, 消去したい節点変位 (slave) をd s として, それに対応して剛性行列と外力も並べ換えを行って K K s K K s d í d s F í Fs ここで従変位に対応する力が解析において無視できるとして F s とおくと 式 (6.b) から [ K ]{ } + [ K ]{ d } { F } s s (6.) d (6.a) [ K ]{ } + [ K ]{ d } { } s d ss s (6.b) { } [ ]{ } d (6.) s K s d となるのでこれを式 (6. a) に代入して ([ K ] [ K ][ K ] [ K ]){ d } { F } となる. よって縮約された自由度の小さい剛性行列は となる. s [ K] [ K ] [ K ] [ K ] [ K ] ss s s ss s (6.5) (6.6) この縮約法は静的な剛性方程式 (.) に基づいているので静的縮約 (statc condensaton) とも呼ばれる. 消去する自由度の選び方が問題であるが, 統一的な基準はなくエンジニアリングセンスで選ぶ. 縮約を使うと解析自由度は小さくなるが, 剛性行列が当初もっていた疎行列性が失われるので計算効率の面で特別有利な方法ではない. この縮約法は振動問題にも適用できてそのときの質量行列は式 (6.) を使って M M M K K (6.7) [ ] [ ] [ ] [ ] [ ] となる. この縮約法を用いた振動解析例を表に示す. s ss s 69

10 静的縮約という言葉があれば当然のこととして動的縮約 (dynac condensaton) なる方法があって衝撃解析などの速い現象に有効であるといわれているが, 一般には使用されない. 6. 減衰 質点系の振動方程式は k u + cu + u f であるが, これと同じく構造においても減衰行列 [C] を導入して (6.8) { } [ K]{ d} + [ C ]{ d} + [ M]{ d } F( t) (6.9) として減衰項を含んだ振動方程式を導入する. 減衰行列は計算によってつくることはできない. 減衰が小さい場合には [ C] a[ M] + b[ K] (6.) とおいて行列の形を [K] か [M] と同じにするとのちのち都合がよい. この形の減衰をレーリー減衰 (Raylegh Dapng) という. 実際には a か b のどちらか一つの項のみ採用し, その値は実験的に決め る. この方法については第 7 章にて述べる. 実際問題として 減衰マトリックスの決定は難しい 構造物においては 粘性減衰よりも乾燥摩擦や履歴減衰 ( 内部摩擦 ) の方が重要となる場合が多いからである 従って通常巴 質量マトリックスの α 比と剛性マトリックスのβ 比を組み合わせて 近似するのが一般的である : [ C] a[ M] + b[ K] 6

11 + c + k f wt f f e {( k w ) + wc} f l C 7 C ± æ Cc ö ç è ø c æ C ö ç è ø k k C c k C a + bk a b + C k k c a b + w w 質量減衰 (α) α 減衰を用いると減衰比は振動数に反比例する つまり 振動数が小さい場合は減衰が大きくなり 振動数が大きい場合は減衰が小さくなる 個の値しか入力できないので α の計算には最も振動数を用いるべきである α の値は 以下のメニューで設定される. 質量マトリックス乗数 [APHAD] として指定する 構造減衰 (β) 構造減衰とは構造物固有の減衰である 多くの問題では 次式のように剛体として粘性減衰は無視される その結果 または ここで w pf f I 固有振動数 a b bw w pf β 減衰を用いると 減衰比は振動数に比例する つまり 振動数が小さい場合は減衰が小さくなり 振動数が大きい場合は減衰が大きくなる あるケースでは 周波数範囲に渡って一定の減衰が要求される場合がある つの減衰関数の合計巴 それらが交差する付近の周波数範囲でほぼ一定となっている 従って 減衰比 (ξ) と周波数範囲 (( 式省略 ) と ( 式省略 )) が与えられれば つの式を同時に満たす α と β を求めることができる k 6 α

12 減衰比 a w + bw a pf a pf + bpf + bpf ( あるケースでは [ ] 6.5 時系列解析 時系列解析とは振動方程式 f f C を用いて一定減衰比を表すこともできる ) u v ku + cv + v f u v ( t) k c u o u í { f ( t) } o í + v o v ku + cu + u f t (6.) ( ) を直接時間軸上で解いていくものである. 加速度は速度の微分, 速度は変位の微分であることを微分を差分に置き換えて計算していくものであり, 差分の取り方により種々の方法がある. 常微分方程式の一般解法として RungeKuttaGll 法があるが, 振動方程式 ( 階の常微分方程式 ) であることを利用して一般的に使われているのは Newark の b 法 Houbolt 法 Wlson の q 法 の つであるが,Newark の b 法 ( しかも b /6 とした線形加速度法 ) が使われることが一番多い. ま ず, 線形加速度法から説明する. ( 全体 ) (β 減衰 ) (α 減衰 ) 6.5. 線形加速度法 まず, 変位を u(t) として微小時間増分を D t とすると u(t+) はテーラー展開により 6

13 u ( t + ) u( t) + u ( t) + u ( t) + u ( t) + ( t + ot) u ( t) u u ot と無限級数となる. これを有限項で打ち切るため u ( t + ) u( t) + u ( t) + u ( t) + u ( t + ) 6 6 (6.a) として加速度の微分に現在の値を使い, その加速度は時刻 tから t+ の間に線形に変化するとして速度は u ( t) + u ( t + ) u ( t + ) u ( t) + (6.b) æ ö ( u + V ) k çu + u + u + V + Cíu + + V f( t) è 6 ø æ ö æ ö ç + C + k V f kçu + u + u è 6 ø è ø æ ö Cç u + u è ø 加速度は振動方程式 (6.) に速度と変位を代入して となる. u ( t + ) ç + c + k [ f( t + ) cíu t 6.5. Newark のb 法 æ è ( ) + u ( t) k u( t) + u ( t) + u ( t) 線形加速度法を一般化すると ot u ( t + of) u ( f) + ( u ( t) + u ( t + ot) ) u t + u t + u t + u t + b u t + ( ) ( ) ( ) ( ) { ( ) u ( t) } í 6 ö ø (6.c) (6.a) となる.b /6 とすれば線形加速度法である. 数値的収束性が良いのは b / のときであり, この / と /6 がよく用いられる. まとめれば u æ ç è ö ø ( t + ) + c + b k f( t + ) c u ( t) + u ( t) æ çè æ ö ( t) + u ( t) + b ( u ( ) kíu ç t è ø (6.b) u ( t) + u ( t + ) u ( t + ) u ( t) + (6.c) í 6

14 である. 6.6 モ ド重畳法 振動問題をその系に特有の固有振動特性から分析することを広義にモード解析という. 固有振動解析はモード解析の中でも最も重要なものである. 振動問題はほとんどがモード解析に終始するといってもよい. その中で, 固有振動モードの重ね合わせで任意の応答を解析する方法をモード重畳法と呼ぶ モーダルパラメータ 振動方程式 [ ]{ d} + [ M ]{ d} { F( t) } K (6.) について考える. 固有振動解析が既に行われて固有振動数 wn と固有振動モード f n がわかっているとき, d t を固有振動モードの重ね合せで 応答 ( ) とする. ここに固有モード行列 [F ] は { ( t) } f q ( t) [ F] { q ( t) } d å l である. 振動方程式に式 (6.5) ここに [ F] [ f, f, f, f ] (6.5) (6.6) n を代入し, さらに左から [ ] F をかけると [ K ~ ]{ q} + [ C ~ ]{ q} + [ M ~ ]{ q } [ F] F [ K ~ ] [ F] [ K] [ F] [ C ~ ] [ F] [ C] [ F] [ ~ M] [ F] [ M] [ F] (6.7) であるが,[K] と [M] については固有ベクトルの直交性からすべて対角行列となり,[C] については [ C ] a[ M] + b[ K] (6.8a) (6.8b) (6.8c) (6.) であるのでやはり対角行列となる. それらの対角成分は k, c, でそれぞれモード剛性, モード減衰, モード質量と呼ばれる. このうちモード剛性とモード質量とは独立でなく固有 ( 角 ) 振動数 w を介して の関係がある. 各対角成分を取りだせば k w (6.9) k q [ t] + c q + q f f (6.) なる n 個のそれぞれ独立した 自由度の振動方程式となる. したがって 6

15 c pw V (6.) となり, モード減衰はモード減衰比 V に置き換えられる. 低次より n 組の 固有振動数 w, 固有振動モード f, モード質量, モード減衰比 V をモーダルパラメータという. モーダルパラメータは外力に無関係な, その構造系が独自にもっている動特性である. また, 固有振動モードはそもそも相対値であり, その絶対値はモード質量と関係する. 式 (6.8) から計算されることからわかるように, k, c も相対的な値であるが, どれかつ決めれば他のつの値も決まる. この固有モード ( ベクトル ) の表示方法として通常次のつの方法が使われている.. 固有振動モードの最大変位をとする.. 固有振動モードのノルムをとする.. モード質量がとなるよう固有振動モードを決める. 番目の方法で決められた固有振動モード f は 番目の方法での と f との間に f f (6.) の関係がある.C 式 (6.) は 自由度系として簡単に解けるのでその線形重ね合わせとして多自由度としての応答が式 (6.5) によって得られる. これがモード重畳法で採用されるモード数は多くの場合 個以下, 多くても 個程度である. なぜなれば振動試験を検討されのものがその程度だからである [6.]. 有限要素ハンドブック, 培風館 参考文献 [6.].Bathe,Wlson( 菊池訳 ) 有限要素法の数値計算米国での有限要素法の振動関係の標準的教科書 [6.]. 戸川隼人 : 有限要素法による振動解析 サイエンス社 (975) 有限要素法そのものでなく, 固有値解析アルゴリズムなどに詳しい [6.].R.R.Crag,Jr.:Structural Dynacs, An Introducton to Coputer Methods, Jhon Wley & Sons (98). 有限要素法の部分構造法に詳しい [6.5]. 小松敬治 : 大規模柔軟構造物のモデリング, 計測と制御 Vol.8,No. (989) pp 図 6. の出典 65