第 章有限要素法 ( その ). 梁要素 有限要素法においては外力も境界条件も節点で考える. もちろん分布荷重は考慮でき るが, 要素上の分布荷重は適当に節点への等価は集中荷重として置き換える. こう考える と梁の曲げの方程式 (.8) において分布荷重無し (p()) の d d (.) である. この一般解は先に解いたように (.) となる. 梁の有限要素の長さを とすると, その両端, にて境界条件を導入して解い てきた. いまここで考えるのは梁の中の一部をきりとった要素である. 未知数が つ ( ~ ) であるので両端で つづつの境界条件を入れればよい. 入れ方として d M or d F or の つの量の組合わせが考えられるが, 右と左の節点で採用する物理量が異なっては統一 性がとれない. 従ってひとつの方法として と d/d の変位を採用する方法があり, こ れを変位法という. さて, この境界条件, d, d d at d at
を入れて,,, について解けば ( ) ( ) 得る. よって となるが ( / ( ) ( ) と無次元化している ), これを { N ( ) } { d} と表現して内部変位 と節点変位 { d} とを結びつける { } { ( ) ( ) } (.) (.) 節点での F と M は図. の符号と図. とを見比べて d d F F F ( ) ( ) F d M d M M ( ) ( ) M N を変位関数と呼ぶ.
であるので, 節点での値を添字, で表すと F M F M となる. これを { F} [ K]{ } と表示して剛性方程式と呼び,[K] を剛性行列という. (.5) (.) 対称行列になっていることに留意されたい. これは偶然でなく, 物理的には MallBtti の相反定理によるものであり,. 節のエネルギー原理による導出によ り明らかになるであろう. 一定の力 F に対して K が大きければ変形 は小さくなるので K が大きいほど剛な構造ということができる. 分布荷重がなければこれが厳密解である. 従 って, 分布荷重のない一様な梁の場合, いくら長くても要素分割は つでよい. 分布荷重 がある場合, 一つの梁をいくつかの要素に分割して分布荷重を等価な集中荷重に置き換え る. これについては次の例題で実例にそって考える... 例 片持梁 左の固定点において で剛性方程式は F, 右の自由端で, M F であるの
M F (.7) である. この式を計算すると上半分と下半分について剛性行列が M F (.7) (.8) となる. 上半分の式 (.8) は変位, が得られたあと反力を計算する式で, 下半分の式 (.9) は外力が与えられたときに変位, を計算する式である. この後者 (.9) について未知数, について解くと ø ö ç è æ となる. また節点 での反力 M, F は得られた節点変位, を剛性方程式の右上半分の式 (.8) に適用して M F となる. これらの結果は前章の例 の厳密解 (.5) と一致する... 例 : 両端単純支持この例では中央に集中荷重がある. 有限要素法では力と変位は節点においてのみ与えられるから中央点にも節点を設け, 要素 節点でモデル化することになる. 節点 と, 節点 と とで作られる剛性行列を重ね合せ, 力の条件 M M M, F と変位の条件 を代入すると
5 図.: 両端単純支持梁 F F F を得る. ここで節点 において両要素が重なっていることに注意されたい. 例 の式 (.9) と同じく未知節点変位に関して第,,, 番目の行と列を抜出して剛性方程式を再構成すると 8
この式を変位について解けばよいが, 文字計算は面倒なので数値例として したがっ ( て /) を採用して計算すると p 8 を得る. 前章の同じ例題の結果 (.8) で とおけば ma 8 となって厳密解に一致する. また となって対称性も満たしている. ここでは表示 上, 剛性方程式を解くのに逆行列を使っているが, 実際には連立 次方程式を解けばよい... 例 : 分布荷重のかかる片持ち梁 p 分布荷重がかかっているときにはモデル化するわけにはいかない. なぜなれば有限要素 を導出するときに採用した変位関数は分布荷重がないとして導かれたものだからである. よってこの例題の場合, 分布荷重を節点上の集中荷重に置き換えてやらなければならない. この例題の場合, つの要素に分割して, 分布荷重の合計 p をその / づつの要素の両
7 端に振り分けて図のように等価な力をかけることにする. として数値計算を行う. 先ほどの つの例と同様に未知変位に関する行と列を残して 5 5 8 sm. 8 8 8 p が得られる. この剛性行列をみればわかるように対角の近傍を除いてゼロである疎行列となっている. 実用的な有限要素法プログラムはこの行列の性質を使って記憶容量と演算回数の節約を計ったコード ( バンド行列法やスカイライン法 ) になっている. この計算結果は.5.75.75 9.5.5.5. 8 p 5 5 である. ここで自由端の変位 5 は.5p/ で, (.7) において とおいて得られる厳密解 p/ に 要素でよく一致していることがわかる.. 立体骨組解析ここまでは簡単な直線梁について考えてきた. 本節では任意の骨組構造の有限要素解析を説明する. 各梁要素が同一直線上になければ曲げ変形は伸び縮みや捩れになって伝わっていく. 立体梁要素の変位とそれに対応する力を図と表とに定義する. 要素の座標系は 軸に沿って梁を置くものとする. 添え字の は節点 の値であることを示す.
方向の変位 X 方向の力 方向の変位 Y 方向の力 方向の変位 Z 方向の力 β 軸回りの捩れ角 M 軸回りのモーメント β 軸回りの回転 M 軸回りのモーメント β 軸回りの回転 M 軸回りのモーメント æ ç è,, ö ø 図.5 立体梁要素 W 次元梁要素を作るためには, 平面における曲げ, 平面における曲げ, 軸方 向の伸び縮み, 軸まわりの捩りを考える. まず, 平面における曲げは式 (.5) より Z M Z M 平面における曲げは Y M Y M (.) (.) 8
下添字の, はそれぞれ 軸, 軸回りの値であることを示す. 次に立体要素として 節 点 自由度を確保するために梁の伸縮要素 ( 棒要素 ) と捩り要素をつくる... 棒要素と捩り要素 棒要素 EA (.)~(.) より F As A( E) 図. EA の解を として ( ) 棒要素 () の条件で EA, を求めると ( ) となる よって F ( ) で F( ), F( ) とすると となる 捩り要素 X X EA GJ 図.7 棒の捩り要素 q 棒の場合と同様に (.) より M GJf GJ C として, を求める操作を行うと GJ (.) 9
GJ M M (.) を得る.. 次元梁要素式 (.)(.)(.)(.) を総合して 次元梁要素の剛性方程式を組み上げると [ ] K M M M Z Y X M M M Z Y X (.5) ここに [K] は EA EA GJ GJ
EA EA GJ GJ (.) と与えられる. これが章末のプログラムの中のサブルーチン fntion[sk]skmat(,ea,z,y,gj) にプログラム化されている さて, 多数の要素で構成される構造の解析においては, 各要素が局所座標系で作られているのでその剛性行列を全体座標系に変換してから例 のような全体剛性行列を作る. 変換行列 [] により局所座標系と全体座標系 ( バー付きで表示 ) の関係は [ ] [ ] [ ] [ ] (.7) ここに [ ] は 掛けるの行列であり
[ ] osj sin j sin j m/q osj n / q m /q mn/q n q (.8) である. ただし X Z 図.8 局所座標と全体座標 (,, ) を節点 の座標,( ) を節点 の座標として ( )/ m ( ) / ( )/ q n m である. 部材の 軸が全体座標系の 軸と平行な場合には である. 式 (.7) を [ ] osj sin j sin j osj { } [ ]{ d} Y d (.8) と書き表せば全体座標系での要素の剛性行列は式 (.) より [ K] [ ] [ K][ ] (.) となる. 章末のプログラムではこの変換行列がサブルーチン fntion[h,fi]hnka (i,j,i,j,i,j) で在る 式 (.) の変換操作はメインプログラムの lmnt
assml 部分で存される ここでの座標変換は実は情報不足で, 実際には I 型梁, 型梁のように 方向曲げと 方 向曲げが異なるような場合の要素の 方向の回転を一意的に表示できない. 系か ら XYZ 系への変換には 軸の方向余弦, すなわち 9 個の情報が必要となるのに i,j 座 標の 個では不足である. これに対応するにはコードアングルの導入 [.] や第 の点を指 定して向きを決める方式がある.. 有限要素法解析の流れ. データ入力 (a) 節点数 () 要素数 () 節点の座標 (d) 要素を構成する節点の構成 () 各要素の,EA,GJ (f) 変位の境界条件 (g) 外力. 各要素での剛性行列の作成. 各要素の全体座標系への変換. 全体剛性行列に組み入れ 5. 変位境界条件の導入. 連立方程式を解く 7. 求められた変位から反力を求める 8. 求められた変位から要素の応力を求める. 一般定式化.. 有限要素法における関係式 変位ー節点関数 { } { N}{ d} (.) 梁の場合, 式 (.) 歪ー変位関係式 { } [ B]{ d} (.)
この関係式は変位節点関係式を微分することにより得られる. 梁の場合 p. で説明した ように d r d { ( ) ( ) ( ) } 応力ー歪関係式.. 仕事 ここに 歪エネルギー とすれば (.) 外部力による仕事 物体力による仕事 梁要素の場合 ò { } [ D] { ε} σ (.) s E である. { } { s} dv { s} [ B] [ D] [ B] dv { d} [ K] [ B] [ D] [ B] d ò ò (.5) { d} [ K] { d} { d} { F} (.) U (.7) ò { } { p} dv { d} { N} { p} dv { d} { F p } U ò (.8) よって物体力と等価な節点力 F は p { F p } ò{ N} { p}dv (.9)
振動などによる慣性力も物体力で, 時間に依存する力は d p dv dt { } とおくことができる. したがって慣性力の場合 ここに である... エネルギー原理 { } [ M] F p { } dm (.) { d} d (.) dt [ M] { N} { N} rdv ò (.) 物体力 p と集中荷重 F を受け, 物体内部で変位 { } と節点変位 { } d とを生じる構造物 の歪エネルギーを とすると, 荷重による外部仕事は歪エネルギーの変化に等しくなけれ ばならないので ( ) d ( U d F ) d (.) p すなわち与えられた変位に対して全ポテンシャルエネルギーが停留値をとる. 歪エネルギ ー は常に正の値をとるのでこの停留値は最小値である. これがポテンシャルエネルギー 最小の原理である. 系の平衡状態では変位に関してしてポテンシャルエネルギーは最小値をとるので, そのと きにはポテンシャルエネルギーを変位 d で微分した値がゼロとなる. すなわち K d F (.) [ ] { } { } { } として剛性方程式が得られる. 動的な運動を考えるときは [ K] { d( )} [ M] d ( t) F p { } F( t) { } t (.5) でこれは振動方程式となっており, 先ほど定義した [M] は質量行列 (mass matri) と呼ば れる... 三角形平面要素 前節の有限要素法の定式化にしたがって演習問題として三角形平面要素を導いてみる. 図.9 三角形平面要素 5
変位 節点関係式 (.a) 5 (.) とする. ここでは核節点で 方向と 方向に 自由度 全 節点で計 個の自由度があるので未定係数を ~ の 個としてある これらは完全 次多項式となっている 点 の (,) 座標を (, ) とし そこの変位を, とすれば なので ( ) { }, であり 同様に ) (, も得られ逆行列を計算して a a a a S (.7) を得る ここに a
7 dt S で S は三角形の面積となっている. a などの添え字はサイクリックに変わっていくものとする. 歪 変位関係式は歪 { } が { } g (.8) として先ほどの変位関数を微分することにより { } [ ] { } d B g S (.9) が得られる. 薄い板が平面内で引張られるような平面応力状態では, 応力 歪関係式は { } [ ] { } s D g n n n E t s s n (.) である. よって剛性行列は [ ] [ ] [ ][ ] [ ] [ ][ ] B B B B K ò D hs dv D (.) として積分なしで得られる.
参考文献 有限要素法の具体的なプログラム作成について [.] 三本木茂夫, 吉村信敏 : 有限要素法による構造解析プログラム 考え方と解説, コンピュータによる構造工学講座 IB, 培風館 97. [.] 小堀為雄, 吉田博 : 有限要素法による構造解析プログラム, 丸善,98. [.] 山田嘉昭監訳 :CA9/SSAN で解く有限要素構造解析 ( フロッピディスク付 ) 丸善, ( 円 ) [.] 信原泰夫, 桜井達美, 吉村信敏 : 有限要素法のプログラムデザイン, コンピュータ による構造工学講座 ⅡB, 培風館 97, pp.8 88. 8