離散的な弾性構造系の静的平衡状態の安定性

4. 構造座屈現象の基礎

4.1. 構造座屈の分岐力学理論

4.1.1. 離散的な弾性構造系の静的平衡状態の安定性

まず，動的な運動状態における Liapunov の安定性を考える．構造物の節点変位ベクトルを時間tの関数としてD( )t { ( )}D t_i ，その速度ベクトルをD( )t { ( )}D t_i とする．これを用いて次の状態空間を定義する．

( ) ( ) ( ) t t

t Z D

D (4.1) そして，自律系

( )t ( ( ))t

Z G Z (4.2) において離散系構造物の運動方程式

( )

MD CU F D P (4.3) を考える．ここでPは時間tに依存しない静的外力ベクトル，Mは質量マトリクス，Cは減衰マトリクス，F D( )は等価節点力ベクトルである．運動方程式(4.3)を状態空間で表すと，次のようになる．

1( )

D D

Z D M CD F P (4.4) なお，節点変位ベクトルDに対応する歪エネルギをU( )D とすると，等価節点力ベクトルは次式で与えられる．

( )

F D D (4.5) 次に，全ポテンシャルエネルギ ( )D を，歪エネルギU( )D と外力による仕事D P^T を用いて

( )D U( )D D P^T (4.6) で定義する．ただし，外力による仕事は保存力であり，外力の大きさや方向が変形形状に依存しない弾性保存系を考える．Liapunov関数の候補として，次のような全エネルギ

W ( ) Z

^{を考える．}

( ) 1 ^T ( )

W Z D MD D (4.7)

なお，Liapunov関数W( )Z は，以下の性質を満たす関数として定義される．

条件1: 状態空間の原点Z 0を含む領域内で，W( )Z およびその導関数は連続である．

条件2: W( )0 0であり，内での原点以外の点でW 0である．

条件3: W( )Z 0である．

式(4.7)の第1項は運動エネルギであり，W( )Z は明らかに条件1を満たす．運動方程式(4.3)で表される系

に対しては，式(4.5)を用いると，次式が成立する．

( )

T T T

T T

W Z D MD D F D P

D MD F P D CD

(4.8)

通常，減衰マトリクスCは正定値行列なので，W( )Z はZ 0のとき負であり，条件3を満たす．Mが正定値行列とすると，式(4.7)の右辺第1項はDの正定値関数である．従って ( )D がDの正定値関数ならば条件2を満たす．更に，次のような剛性マトリクスK を定義する．

2 2

i j

U U K F

D D (4.9) 運動方程式(4.3)で表される系の静的平衡状態の特性について定常状態の安定性から考える．いま，自律

系(4.2)の一つの状態

Z Z

₀で

( )0

G Z 0 (4.10) が成立するものとする．このとき，式(4.2)より

Z 0

^なので，

D D 0

^であり，

D

^{は時間の経過に}

伴って変化しない．また，

t

^{に依存しない外力}

P

の作用する運動方程式系では，式(4.4)より，

Z Z

₀において

F P

が成立する．従ってこのような状態は静的な平衡状態である．安定性の検定のために微小な振動を加えることを考えると，

F D ( ) P

^，および

P

^が

D

に依存しないことを用いて，運動方程式の状態空間表現式(4.4)の右辺をD 0^において

D

^{に関して線形化し，式}(4.9)を用いると，

1 1

1( )

D 0 I

Z DD M CD F D M K M C Z

(4.11)

となる．式(4.11)は自由振動の方程式の状態空間表現である．右辺の行列を

1 1

0 I

S M K M C (4.12) と定義する．Sは非対称行列であるので，その固有値は一般に複素数である．S^{の全ての固有値の実部} が負である時，原点Z 0は漸近安定であり，自由振動は原点に収束する．いま，

M

^とC^{が正定値と} すると，

K

が正定値であれば上記の漸近安定性に関する条件は満たされる^[39]．

4.1.2. トラスアーチでの構造座屈の計算例

2本のトラス要素で構成されるFig. 4-1のようなトラスアーチを考える．式(4.6)に示した通り，全ポ

( , )D f (4.13)

例えば，Fig. 4-1に示すトラスアーチでは，節点③に赤矢印で示すような静的外力Pを加えたとして，

EAf

P e

(4.14) と書ける．ここで，Eはヤング率，Aはトラスの断面積，e_xは x 軸方向の単位ベクトルである．また，

節点変位ベクトルDとして，は3つの節点のx y z, , 方向の変位を縦に並べたベクトルとして

1 2 9

[D D ... D ]^T

D (4.15) と書ける．従ってFig. 4-1に示すトラスアーチにおける全ポテンシャルエネルギは

1 9

( , )D f ( , ... ,D D f, ) (4.16) と書ける．式(4.13)について，全ポテンシャルエネルギ停留条件は次式で得られる．

( , )

F D 0

D (4.17)

Fig. 4-1に示すトラスアーチでは，

1 9

( , ... , , )

( , ... , , ) 0 , 1,2,...,9

D D f

D D f with i

F D (4.18)

と書ける．そしてこれが解くべき離散化された平衡方程式，即ち静解析における剛性方程式である．

ここで，式(4.17)の釣合点での安定性を調べる．ある釣合点( , )D f の近傍における，f 一定の下での，増分変位 Dに対する全ポテンシャルエネルギの増分量 ( , )D f は次式で得られる．

( , ) ( , ) ( , )

1 2 1

f f f

D D D D

D D F D

D D

D K D

(4.19)

ただし，釣合点での釣合条件式(4.17)と剛性マトリクスの定義式(4.9)を用いた．Kの固有値を _i，固有ベクトルを _iと表す．ただしi 1, 2,...,nであり，Fig. 4-1に示すトラスアーチでn 9である．ポテンシャルの存在により保証されるKの対称性により，この固有ベクトル _iが正規直交系をなすように選ぶことができる．これを用いて直交行列

1 2

[ ... _n]

H (4.20) を定義し，この行列による座標変換

1 2

, [ ... ]

n T

i i n

v with v v v

D H v v

(4.21)

を用いて式(4.19)を書き直すと，次式を得る．

2 1

( , ) 1 2 1 2 1 2

T T

n n

i i i

D v H KH v

v v (4.22)

この式から， ( , )D f が極小値を取るための十分条件は，すべての固有値 _iが正となることであると分

かる．即ちこの系の安定性は

K

の固有値が全て正安定 (4.23)

K

の固有値が一つでも負不安定 (4.24) と判定される．これは4.1.1節の安定性の論旨と一致する．

それでは実際にFig. 4-1に示すトラスアーチにおいて具体的に計算してみる．行列で記述すると煩雑になるため，拘束条件を加味し，スカラーの式で記述する．即ちこの例では節点③に

x

^{方向の荷重を載}

荷し，節点③の

y z ,

方向は拘束する．節点①と節点②は完全拘束されているため変位が生じない．よってこの系は 1 自由度であるので，節点③の

x

^{軸方向の変位を}d^{とおき，変数を}d^と

f

^{で記述していく．}

また，Fig. 4-1における

H

^と

W

^は，

H 0.05[ ] m

W 1 [ ] m

^とし，

EA 2.8 10 [N]

⁷ とした．

変形前の要素長は W² H²，dだけ変位した時の要素長は _W² ₍_d _H₎²であるため，全ポテンシャルエネルギは

2 2 2 2 2

2 2

2 2 2 2 2

2 2

( , ) ( )

1 ( )

d f EA W d H W H EAfd

W H

EA W d H W H fd

W H

(4.25)

式(4.25)をd に関して偏微分することで，釣合式が次式のように得られる．

2 2 2 2

( , )

1 1

2 ( )

( )

0 F d f

EA d H f

W H W d H

(4.26)

式(4.26)を解くと，荷重変位曲線は次式のように得られる．

2 2 2 2

1 1

2 ( )

( )

f d H

W H W d H (4.27)

更に，式(4.26)をdに関して偏微分することで，接線剛性が次式のように得られる．

2 2 2 2 3/2

( , ) 2 1

( )

K d f F d EA W

W H W d H

(4.28)

座屈点は，式(4.27)と座屈条件式K 0^より，

H 0.05[ ] m

W 1 [ ] m

のとき次のように得られる．

1/3 1/2 1/3 3/2

(d_cr, f_cr) 1 (2 1) , 2(2 1) (4.29) 荷重-変位曲線はFig. 4-3に，固有値-変位曲線はFig. 4-4に示す通り得られる．固有値

K

^は

F

^{の変位に対}

解はFig. 4-3の赤破線矢印で示すように飛び移る．なお，本モデルにおいて構造座屈は部材座屈より小さい荷重で生じる．部材座屈を起こすオイラーの座屈荷重がおよそ9188[N]であるのに対し，構造座屈が生じるときの載荷荷重はおよそ1344[N]となっている．このことから，構造物の座屈は必ずしも部材座屈だけに注意していれば良いというわけではないことが言える．

最後に，分岐座屈について述べる．分岐座屈が起きるのはFig. 4-2に示すような，節点③の

y

^方向が

拘束されていないトラスアーチモデルにおいて，

H W

となるような細長いトラスアーチの場合である．これは 2 自由度の問題であり，式がやや煩雑になるため，また，後述する通り静解析と同じ結果が得られるため，ここでの具体的な計算は行わず，次節にて示す．

Fig. 4-1 トラスアーチのモデル

Fig. 4-2 分岐座屈を生じるトラスアーチのモデル

W

H

①

② ③

x y z

EAf

要素1

要素 2

H 0.05 [m]

W 1 [m]

EA 2.8×10

⁷

[N]

拘束条件

節点①のx,y,z 節点②のx,y,z 節点③のy,z 自由度数 1

飛び移り座屈

W

H

①

② ③

x y z

EAf

要素1

要素 2

H 3 [m]

W 0.01 [m]

EA 2.8×10

⁷

[N]

拘束条件

節点①のx,y,z 節点②のx,y,z 節点③のz 自由度数 2

分岐座屈

Fig. 4-3 分岐理論により得られる飛び移り座屈の荷重-変位曲線

Fig. 4-4 分岐理論により得られる飛び移り座屈の固有値-変位曲線

-1.0 -0.5 0.0 0.5 1.0 1.5 2.0

0.00 0.02 0.04 0.06 0.08 0.10 0.12

分岐パラメータ f[ × 10 -4 ]

変位量 d[m]

不安定

-10 -5 0 5 10 15 20 25 30 35 40

0 0.02 0.04 0.06 0.08 0.1 0.12

固有値 [ × 10 4 ]

変位量 d[m]

4.2. 構造座屈の有限要素法による構造解析

Fig. 4-1に示したトラスアーチのモデルに対して，荷重制御法による静解析を行った結果をFig. 4-5

と Fig. 4-6 と Fig. 4-7 に示す．静解析では式(2.66)，(2.67)，(2.63)を用いて Newton 法にて解を得る．

0.05[ ],

H m W 1.0[ ], m EA 2.8 10 [N]

⁷ ^{とし，ステップごとに}

f

^を

2.5 10

⁸^{ずつ増加させながら}

載荷し続ける条件で解いた．また，応答図は見やすくするためにx軸方向のスケールを20倍に拡大して表示している．Fig. 4-5より分かる通り，Step1919からStep1920の1ステップでいきなり座屈後形状に飛び移っている．前節でも述べたとおりであるが，Fig. 4-7を見ても分かるように，静解析では不安定領域の解の経路を知ることができない．それは，Newton 法の反復計算で _xを更新する際に，

1( )

x x K P f としているため，剛接マトリクスKの固有値に一つでも0が表れるとKは非正則となり，計算が止まってしまう．しかし，数値計算においては離散的なステップで釣合式を解いていくことや，桁落ち誤差が生じることから，厳密に0という固有値が出ることはまれである．本解析例でもFig.

4-5，Fig. 4-7から分かる通り，微小な正の固有値を得た次のステップで突然正の大きな固有値に飛んでい

る．もし微小な正の固有値が出た次のステップで分岐パラメータを減少させたとしても，不安定経路に行くのではなく安定領域側に経路を戻っていくような解経路をたどることになる．よって荷重制御静解析では負の固有値を得ることはできず，不安定領域の解は得られない．但し，弧長制御静解析を行えば，

不安定領域での解の追尾は可能であることが知られている．

次に，同じトラスアーチのモデルに対して動解析を行った結果をFig. 4-8とFig. 4-9とFig. 4-10に示す．動解析では式(2.72)，(2.73)，(2.67)，(2.63)を用いて Newton 法にて解を得る．静解析と同じく

0.05[ ],

H m W 1.0[ ], m EA 2.8 10 [N]

⁷ ^{とし，ステップごとに}

f

^を

2.5 10

⁸^{ずつ増加させながら}

載荷し続ける条件で解いた．また，時間ステップ幅はdt 5 10 [ ]⁵ s とし，速度および加速度の初期値は 0とした．静解析同様，応答図は見やすくするためにx軸方向のスケールを20倍に拡大して表示してい

る．Fig. 4-10を見ると，静解析に見られたような 1ステップでの形状変化は見られず，飛び移り座屈の

途中の形状も得られている．これは，式(2.60)および式(2.67)，(2.73)を比べることで分かるが，動解析のヤコビアンには慣性項が入っており，慣性項に表れる質量マトリクスが正定値の正則な行列であるため，

ヤコビアンの正則性が保証されることによる．なお，Fig. 4-10に示す固有値とは，ヤコビアンの弾性項，

即ち剛性マトリクスの固有値である．Fig. 4-9，Fig. 4-10とも，静解析よりも連続的なグラフが得られ，

Fig. 4-10はFig. 4-4に示した分岐理論により得られるグラフと同じ形状が得られた．Fig. 4-9がFig. 4-3と

異なる形状となったのは，載荷条件を時間ステップあたり

2.5 10

⁸ずつ増加させ続けたためである．従

って，Fig. 4-3の赤破線矢印に示した飛び移りの現象が実際に起こっていることがFig. 4-9から確認でき

る．そしてFig. 4-10から分かるように，動解析では不安定領域の解は，負の固有値を検出できたことによって認知できる．

最後に，前節で触れた分岐座屈について説明する．Fig. 4-2 に示したトラスアーチのモデルで，

3 [ ]

H m

^，

W 0.01 [ ] m

^，

EA 2.8 10 [N]

⁷ ^{とし，ステップごとに}

f

^を

2.5 10

⁸^{ずつ増加させなが}

らx軸方向に載荷し続ける条件の下，解析を行った結果を記す．2自由度系であるため，固有値は2つ得られる．また，この時の部材座屈荷重はおよそ1023[N]であるのに対し，分岐座屈を生じるときの載荷荷重はおよそ623[N]となっている．

まず，y 軸方向には一切の初期不正も載荷も無い対称載荷条件にて静解析を行った結果を Fig. 4-11

からFig. 4-15に示す．応答形状は見やすくするためにy軸方向のスケールを20倍にして表示している．

一つ目の固有値Fig. 4-14を見ると負の固有値が表れており，変位3.33×10^-5[m]付近で固有値が正から負に転じている．しかしFig. 4-11，Fig. 4-12，Fig. 4-13にはその点での際立った変化は見られず，応答は安定領域から大きな変化は無く連続しており，x軸方向のみの値の推移が見られる．これは，応答がFig. 4-15 に示すもう一つの固有値解経路を辿っているためである．即ち，Fig. 4-15に示す2つ目の固有値に対応する固有モードは節点③の変位について

( , , ) ( 1,0,0) x y z

^{であるが，}Fig. 4-14に示す一つ目の固有値については，

( , , ) (0,1,0) x y z

となっている．そして，完全対称形である本モデルでは応答が辿る経路は

( , , ) ( 1,0,0) x y z

のモードであるが，もう一つのモード

( , , ) (0,1,0) x y z

は固有値が正から負へ転じるモードであり，固有値が負である不安定領域では微小な外乱や初期不正がある非対称モデルであれば短時間で大きく動くような座屈を起こす状態となっている．また，Fig. 4-13 が単調増加であることから，

一つの荷重に対して複数の釣合点が存在することが無いため，飛び移りのような現象は起きず，これが分岐理論と静解析で同様の結果が得られる理由である．更に，分岐座屈では静解析でも負の固有値を検出することができ，これは即ち動解析でもヤコビアンの弾性項において同様の固有値解経路が得られることを意味する．

そこで次に，節点③のy 軸方向にも

EAf

^の1/80000 という微小な荷重を載荷する非対称載荷条件にて構造動解析を行った結果をFig. 4-16からFig. 4-22に示す．応答形状は見やすくするためにy軸方向のスケールを20倍にして表示している．また，Fig. 4-19とFig. 4-21はそれぞれFig. 4-18とFig. 4-20の

Step1500までの結果を表したものである．一つ目の固有値Fig. 4-21を見ると，対称載荷の静解析とほぼ

同じく，x軸方向の変位3.34×10^-5[m]付近で固有値が正から負に転じている．そしてその後，非対称載荷

はFig. 4-17に見られるように，節点③の急激な変位をもたらす．応答形状も対称載荷とは異なり，Fig. 4-16

のようにy軸方向に倒れる応答を示す．座標に急激な変化が現れる辺りでFig. 4-19に示す荷重-変位曲線にも急激な変化が現れている．非対称載荷では負の固有値が現れるモード

( , , ) (0,1,0) x y z

^{が応答に出現}

しており，微小な載荷荷重の偏りにも関わらずドラスティックな変位を生じさせる．これが分岐座屈と呼ばれる現象である．Fig. 4-22 に示すもう一つの固有値は対称載荷と同様の値を示した．最後に，同様の非対称載荷を静解析で解いた応答形状をFig. 4-23に示す．静解析ではStep889から負の固有値に転じたが，その 1ステップで節点③のy座標が正から負に飛んでしまっていることが確認できる．静解析では慣性項が無いため，部材の内力と載荷荷重が釣り合う形状を探すことになる．従って実際の分岐座屈現象とは異なる不連続な応答形状推移を示す．このことから，静解析では負の固有値を得られることはあるが，固有値が負の領域に入ると実際の現象に見合う妥当な応答を得ることはできないと言える．

なお，動解析において不安定性を論じていく中で，剛体運動を含めた系全体の運動が発散するという動的安定性における不安定と，座屈によってヤコビアンの弾性項の固有値が負になるという構造形状不安定を混同しないようにする必要がある．ここまでの座屈理論では式(4.10)に示す通り，動的安定性に関する概念を除去したアドホックな系で論じている．本研究では構造形状不安定における運動を取り扱っていくため，以降の章でも式(4.10)を前提とした系の中で理論を記述していく．また，形状不安定のことを不安定と呼ぶこととする．

ドキュメント内目次 (ページ 36-55)

4. 構造座屈現象の基礎

4.1. 構造座屈の分岐力学理論

4.1.1. 離散的な弾性構造系の静的平衡状態の安定性

W ( ) Z

Z Z

Z 0

D D 0

D

t

P

Z Z

F P

F D ( ) P

P

D

D

M

K

4.1.2. トラスアーチでの構造座屈の計算例

EAf

P e

, [ ... ]

v with v v v

D H v v

K

K

x

y z ,

x

f

H

W

H 0.05[ ] m

W 1 [ ] m

EA 2.8 10 [N]

H 0.05[ ] m

W 1 [ ] m

K

F

y

H W

W

H

①

② ③

x y z

EAf

要素1

要素 2

H 0.05 [m]

W 1 [m]

EA 2.8×10

[N]

拘束条件

節点①のx,y,z 節点②のx,y,z 節点③のy,z 自由度数 1

飛び移り座屈

W

H

①

② ③

x y z

EAf

要素1

要素 2

H 3 [m]

W 0.01 [m]

EA 2.8×10

[N]

拘束条件

節点①のx,y,z 節点②のx,y,z 節点③のz 自由度数 2

分岐座屈

-1.0 -0.5 0.0 0.5 1.0 1.5 2.0

0.00 0.02 0.04 0.06 0.08 0.10 0.12

分岐 パラ メー タ f[ × 10 -4 ]

変位量 d[m]

不安定

-10 -5 0 5 10 15 20 25 30 35 40

0 0.02 0.04 0.06 0.08 0.1 0.12

固有値 [ × 10 4 ]

分岐パラメータ f[ × 10 -4 ]