骨組構造物の非線形挙動の数値解析法に関する研究
建設構造研究室 松野 純一
指 導 教 官 岩崎 英治
長井 正嗣

1

序論

幾何学的非線形解析手法に関する限り，理論的な基礎はほぼ完成していると考えられている. その一方で，実際に非線形解析を必要とする設計者や技術者にとって，複雑な理論や，煩雑なアルゴリズムよりなる数値解析法を理解するのは，時間的制約もあり困難となっている. このため，収束性や求解に必要な，物理的に不明確なパラメータを減らし，入力データが容易に作成できるような非線形解析システムの開発が必要である.
すなわち，弧長の設定に関しては，一般には，釣り合い曲線の曲率の大きな部分では，弧長を小さくし，逆に計算効率の観点から，曲率の小さな部分では，弧長を大きくするべきである．しかし，事前に適切な弧長を入力データにより設定するのは困難であり，試行錯誤によりデータ作成と解析を繰り返しながら適切な弧長を設定することになる．
本研究では，入力データが容易に作成でき，容易に非線形解析を行えるような解析システムを開発するための基礎的研究に相当するもので，汎用性の高い弧長制御法に，不平衡ベクトル最小法 1)を用いて，その増分量である弧長の長さを曲率の変化に伴い自動的に決定できるような新しい弧長自動設定法を提案する．

22.1

数値解析法
非線形平衡方程式

変位法によるマトリクス構造解析法に現れる連立非線形方程式は，構造全体の力の釣り合い条件を表す平衡方程式である．この方程式は次式のようになる．
F (D) = P
(1)
ここに, F は内力ベクトル, P は外力ベクトルで,保存力であるとする．なお，内力は初期状態から現時点までの変位の履歴で表されるが，表現を簡略化して変位 D で表す．また，解析を通して荷重のモードが一定である場合には，ある一定の荷重ベクトル (基準ベクトル)P と荷重の倍率 λ の積で外力ベクトルは表される．このような問題では上式は次のようになる．

F (D) = λP(2)

上式のような n 個の変位成分と荷重倍率 λ で表される平衡方程式の解は，(n+1) 次元空間での曲線を表している．すなわち，上式を解くということは曲線を求めることに他ならない．この曲線を釣り合い曲線と呼んでいる．
上式を増分形で表すと，次のようになる．

K(D)∆D = ∆λP + R(D)(3)

荷重倍率の増分 ∆λ を与えて解 ∆D を求める方法
を増分法，特に荷重を制御変数とする方法を荷重増分法と呼んでいる．なお，上式の R は次式で定義されるベクトルである．

R(D) = λP − F (D)(4)

これは，増分量に関して線形化されたために生じる誤差で，不平衡力と呼ばれる．増分計算により平衡状態の解が得られている場合にはゼロとなる量である．

2.2

弧長制御法

弧長制御法では，変位ベクトル ∆D と荷重倍率∆λ から求められる釣り合い曲線の弧長 ∆S を与えて，変位と荷重を求める．このための方法として，種々の方法が考えられるが，ここでは係数行列の対称性を損なわない方法を述べる．
弧長は ∆S 2 = ∆D T ∆D + ∆λ2 により与えられるが，変位，回転と荷重倍率の次元や数値的な大きさが異なるので，スケールパラメータ αi (i =1, 2, . . . , n) と αp を導入し，弧長を次のように表す．

∆S 2 = ∆D T A2 ∆D + αp2 ∆λ2(5)∆D1ii

∆λacn

∆λ

∆D2i

∆S

∆Sd

∆λ1i
∆λ2i

1n−1K1i

Loadbn−2Displacement

図–2 不平衡ベクトル最小法
図–1 弧長制御法
反復計算中の修正量 ∆L は弧長と同様の式で表される．

ここに，行列 A は次式で与えられる．

A=



α1
α2..

.

2T 222∆L = ∆D A ∆D + αp ∆λ(6)∆D を分解すると，荷重倍率の 2 次式となり，∆Lを最小とする荷重倍率が，微分により次のように得られる．これが，不平衡ベクトル最小法である．

αn式 (3) の解 ∆D を，基準荷重 P による変位 D 0と不平衡力 R による変位 D 1 に分けて表し，これを式 (5) に代入すると次式が得られる．

∆S 2 = (∆λD 0 + D 1 )T A2 (∆λD 0 + D 1 )+αp2 ∆λ2(7)これを荷重倍率に関して解くと，√−Y± Y 2 −X ZX

T 22
 X = αp + D 0 A D 0∆λ =(8)

Y = D 0 T A2 D 1Z = ∆S 2 − D 1 T A2 D 1上式より，弧長の式を満足する荷重倍率は二つあることが分かる．これは，図 1 のように (n + 1)次元空間内での球面と接線との交点が二つあるためである．このとき，荷重倍率としてどちらを用いるかが問題になるが，既に求められた平衡点から離れたほうを用いると良い．
このようにして得た解は，釣り合い方程式を線形化しているために，正確な平衡状態にない．そこで，不平衡力を許容値以下に収める反復解法を用いる必要がある．この手法には既往の研究で幾つか考案されているが，本研究においては，文献
1) でその有効性が確認されている，不平衡ベクトル最小法を用いる．

T∆λ = −3f
D A2 DT

αp 2 + D A2 D

(9)

弧長増分量の自動設定

弧長制御法においては，反復計算により誤差を低減する一方で，さらに，非線形性が強い場合には, 増分量を小さくし, 非線形性が弱い場合には,増分量を大きくすることにより, 反復回数の低減や計算コストの節約を図ることが期待できる．しかし, 釣り合い曲線の形状が事前に分からなければ, 適切な増分量を設定することは困難なので, 釣り合い曲線の非線形性の強弱, 言い換えれば, 釣り合い曲線の曲率に応じて増分量の大きさを, 自動的に設定できるような手法が必要になる.
次に，既往の研究で考案されている弧長増分の大きさを自動的に変更する計算方法を 3 手法に限って紹介する．また更に，本研究にて新しく提案する方法についても説明する．

3.1

Itoh と Nogami による弧長自動設定法

Itoh と Nogami の方法 2) では，求めようとする次の増分ステップにおける弧長の大きさを次式により求めている．
∆sn+1 =

α + Ns Nc∆sn
α + Nc2(10)ここで, Ns は目標収束回数, Nc は前回の釣り合い点を求めるのに要した収束回数，α は任意の定数であり，これは，α = 0 のときには，次の弧長はNs /Nc 倍になり，α を大きくすると次の弧長は，直前の弧長に近い値に漸近するような定数である．
Itoh と Nogami の方法では，弧長の自動計算を始めるために，直前の増分段階での弧長と収束回数が必要なため，最初の増分ステップの弧長 ∆s1と，この他に目標収束回数 Ns と任意定数 α の合計 3 個のパラメータが必要である．

3.2θn∗

荷重∗θn−1nn−1
snn−2

sn−1

図–3 W. Szyszkowski らの方法

Chaisomphob ら 3) は，弧長を釣り合い曲線のC∆sn+1 =
κ(s)

θn
荷重

θ2(11)

vu µ 2¶µ 2¶nud λ(s) 2 Xd Di (s) 22t2κ(s) = α+αds2ii=1Szyszkowski ら 4) は直前に求めた 3 つの解と, 次の二つの仮定• 各釣り合い曲線上の点を通る曲線の曲率は線形変化する• 各釣り合い点で、直前に求めた釣り合い点からこの点を通り，次の釣り合い点まで延長した直線と，次の釣り合い点とのずれは等しいを用いて，次の弧長を次式により求めている.
θn∗

∆sn ,
∗θn+1

∗
∗
θn+1
= 2θn∗ −θn−1sn+1n−1 sθnn θ κn = snn−1
κn−1 = sn−1

κ1 s22 = κ2qs23 = · · · = κn s2n+1
sn+1 = s2 κκn10

Szyszkowski らによる弧長自動設定法∆sn+1 =s1ds2(12)ここで，λ は初期状態からの荷重倍率，Di は変位成分である．上式に含まれる (d2 Di /ds2 ) と (d2 λ/ds2 )を得るために，釣り合い曲線上の複数の点での値から 3 次式で近似している.
Chaisomphob らの方法では，弧長の自動計算を始めるために，直前の 3 点の解が必要なために，最初の 3 ステップ分の弧長 ∆s1 , ∆s2 , ∆s3 と，前述の任意係数 C の合計 4 個のパラメータが必要である．
3.3

κn−1 snn

κ12s2 κ = θ22s2s21 κ1 = θs11θ1

ここで, C は任意の定数である. また，κ(s) は釣り合い曲線の曲率であり，これは次式により求められる．
psn+1

変位

Chaisomphob らによる弧長自動設定法

曲率の逆数に比例するように決めることを提案し，次式を示した．

∗θn+1∆sn∆sn−1

変位
図–4 弧長の自動設定ここで, θn∗ は, 平衡点 (n − 2) から点 (n − 1) を通る線と点 (n − 1) から点 n を通る線の間の角度である.
Szyszykowski らの方法では，弧長の自動計算を始めるために，直前の 3 点の解が必要なために，最初の 3 ステップ分の弧長 ∆s1 , ∆s2 , ∆s3 の合計3 個のパラメータが必要である．

3.4

本研究で提案する弧長自動設定法

図 4 をもとに増分量の設定方法を示す.
増分計算により点 n までの解が得られているものとする. 点 n と一つ前の点 (n−1) の距離を ∆sn ,同様にその前の点までの距離を ∆sn−1 と表す. また, 点 (n − 1) の接線とさらに一つ前の点 (n − 2)の接線の間の角 θn−1 を, これらの点の距離 ∆sn−1で割ったものを点 (n − 1) の曲率 κn−1 とする. すなわち,
θn−1(13)
κn−1 =∆sn−1

同様に, 点 (n − 1) と次の点 n の間の角度変化 θnと距離 sn から曲率 κn は，次のようになる．

κn =κn (∆sn+1 )2 = κn−1 (∆sn )2 = · · · = κ1 (∆s2 )2
これより,∆sn+1 = ∆sn

κn − 1= · · · = ∆s2
κnr

κ1
κn

全ての増分段階での平衡点からのずれが等しくなるように弧長を決めるので, ∆sn+1 として用いる式は, 上式のどれを用いても同じである. しかし, 数値計算上の誤差の累積や, 各増分段階での繰り返し計算の手法として, 不平衡ベクトル最小法などの最初に設定した球面上から外れた点に収束する手法を用いた場合には, 各増分段階でのずれ変位 κn ∆s2n+1 が等しくなっていないので, ここでは, 常に, 最初のずれ変位 κ1 (∆s2 )2 に等しくなるように決めることにする. すなわち, 次式により弧長 ∆sn+1 を決める.
r∆sn+1 = ∆s2

κ1κn(14)

上式に含まれる κ1 , ∆s2 は, 初めの二つの平衡点1 と 2 での解が, 明らかになっていなければ求められない. したがって, ここに提案する手法は, 初めの 2 点については, 弧長を事前に設定する必要があることに注意する必要がある.

44.1

yθn∆sn

この曲率に次の点 (n + 1) までの距離 ∆sn+1 を乗じたもの κn ∆sn+1 は, 点 (n + 1) での角度変化
θn+1 を近似している. また, この角度変化に距離∆sn+1 を乗じたものは, 点 (n + 1) での釣り合い曲線上の点からのずれ変位 (誤差) を表す指標と考えられる. この誤差の指標が常に同じ大きさになるように弧長 ∆s を決めると, 各増分計算での反復回数の均等化が図られるものと予想される.
そこで, この方法を用いて次の点までの距離を決めることにする. すなわち,

sy数値計算例及び考察
トラスドームの有限変位解析

図 5 のようなトラスドームの中心の節点に鉛直
荷重 P ，外周の節点に 2P が作用したときの荷重変x4m

2P2Pz P2P2P

8m

1m
1.3m

E=200GN/m2 , A=1.0m2
図–5 トラスドーム位曲線を考える．このトラス部材のヤング係数はE=200GPa，断面積は A=1m2 としている．弧長増分法の制御変数には，荷重倍率 λ と荷重作用点の鉛直変位 w を用い，これらのスケールパラーメー
タはそれぞれ αp = 0.2, αw = 1.0 を用いる．また，基準荷重は P0 =100MN．収束判定には，収束計算途中の変位の修正量のノルムが各増分段階の最初の変位増分ノルムの 0.001 倍以下としている．
図 6(a) に本研究で提案した弧長自動設定法の計算結果を示す．ここでは，最初の 2 ステップ分の弧長として，∆s = 0.15 を用いている．図示のように曲率の大きな個所では弧長が小さくなり，曲率の小さな個所では大きな弧長となっていることが分かる．図 (b) は，Itoh らによる弧長自動設定法による結果であり，Ns = 3, α = 3.0，最初の増分ステップの弧長は ∆s = 0.15 としている．最初の極大点以降は，弧長が小さくなり，曲率の小さな部分でも弧長は小さなままになっている．図 (c)は，Chaisombphob らによる弧長自動設定法による結果であり，C = 0.2, 最初の 3 ステップ分の弧長は ∆s = 0.15 としている．曲率の大きな部分では弧長は小さくなり，曲率の小さな部分では大きくなっているが，その弧長が非常に大きく変化している．釣り合い曲線の形状が更に複雑な形状となる問題では，弧長が大きすぎると収束解が得られない場合があるので，任意定数の決め方に注意する必要がある．図 (d) は，Szyszkowski らによる弧長自動設定法による結果である．最初の 3 ス
テップ分の弧長は，他の結果と同様に，∆s = 0.15を用いているが，途中から釣り合い曲線を逆戻り

Load factor ( lambda )Load factor ( lambda )2
10-1-22
10-1-20 -0.5 -1 -1.5 -2 -2.5 -3
Displacement ( w )0 -0.5 -1 -1.5 -2 -2.5 -3
Displacement ( w )(b) Itoh と Nogami 法

2Load factor ( lambda )Load factor ( lambda )(a) 提案法10-1-20 -0.5 -1 -1.5 -2 -2.5 -3
Displacement ( w )210-1-20 -0.5 -1 -1.5 -2 -2.5 -3
Displacement ( w )(c) Chaisomphob 法(d) Szyszkowski 法

図–6 トラスドームの計算結果

y1.5mP6mx50m
50m
2
E=200GN/m , A=1.0m2 , I=1.0m4
図–7 扁平アーチ
しだし，釣り合い経路を求めることができなかった．最初の弧長の大きさを変えて計算したが，同じような結果になり，釣り合い曲線を求められなかった．

4.2扁平アーチの有限変位解析

図 7 のような扁平なアーチの頂点から 1.5m ずれたところに鉛直荷重が作用したときの荷重とアーチクラウンの鉛直変位の関係を求める．このアーチのヤング係数は E=200GPa，断面積と断面 2 次モーメントはそれぞれ A=1.0m2 , I=1.0m4 である．弧長増分法の制御変数には，荷重倍率 λ と鉛直変位 v を用い，これらのスケーリングパラメー

タはそれぞれ αp = 0.1, αv = 0.25 を用いる．基準荷重は P =100MN．を用いる．収束判定には，トラスドームの計算例と同様に，収束計算途中の変位の修正量のノルムが各増分段階の最初の変位増分ノルムの 0.001 倍以下としている．
図 8(a) に本研究で提案する弧長自動設定法による結果を示す．これら最初の 2 ステップの弧長を∆s = 0.045 としたときの結果である．図より，荷重変位曲線の曲率の大きな部分では，小さな弧長になり，曲率の小さな部分では大きな弧長となっている事がわかる．図 (b) に Itoh らによる結果を示す．最初のステップの弧長は ∆s = 0.045，他のパラメータは Ns = 3, α = 3.0 を用いている．図より，曲率の大きな箇所を通ると弧長が小さくなり，曲率の小さな箇所に来ても弧長が大きくなっていないことが分かる．これは，弧長が小さくなりすぎたので，途中で計算を打ち切っている．図(c) は Chaisombphob らによる結果であり，最初のステップの弧長は ∆s = 0.045，他のパラメー
タは C = 0.01 とした．図より，λ = 0 付近で，弧長が大きくなりすぎたため，途中で解が発散しLoad factor ( lambda )Load factor ( lambda )

32.5
2
1.5
10.5
0-0.5-1-1.5032.5
21.5
10.50-0.5-1-1.5-2 -4 -6 -8 -10 -12 -14
Displacement ( v )0(b) Itoh と Nogami 法

32.5
2
1.5
1
0.5
0-0.5
-1-1.5

Load factor ( lambda )Load factor ( lambda )(a) 提案法

0-2 -4 -6 -8 -10 -12 -14
Displacement ( v )

3
2.5
2
1.5
1
0.5
0-0.5
-1-1.5-2 -4 -6 -8 -10 -12 -14
Displacement ( v )0(c) Chaisomphob 法-2 -4 -6 -8 -10 -12 -14
Displacement ( v )(d) Szyszkowski 法

図–8 扁平アーチの計算結果た．図 (d) は，Szyszkowski らによる結果で，最初のステップの弧長を ∆s = 0.045 とした．これもChaisombphob らの結果と同様，λ = 0 付近で，弧長が大きくなりすぎたため，途中で解が発散した．

5結論

弧長制御法に基づく不平衡ベクトル最小法に対して， その増分量である弧長の長さを曲率の変化に伴い自動的に決定できるような手法を提案し，既往の研究で考案されている方法との違いを数値計算により比較検討した結果，次のような結論を得た．
(1) 既往の研究で考案された方法では，任意定数等の物理的に不明確なパラメータを必要とするが，その設定方法が不明であったり, そのパラメータにより収束性が大きく変化するために，試行錯誤が必要であることが確認できた．
(2) また, 釣り合い曲線の非線形性状によって, 極端に弧長が短くなったり長くなったりすることがあり, その使用に注意する必要のあることが分った．

(3) これに対して，本研究で提案した方法では，既往の手法と同様に, 最初の数増分ステップ分の弧長を与えるだけで，複雑な形状の釣り合い曲線を持つ構造においても，安定した数値解の得られることが分った．
したがって，本研究で提案した弧長自動設定法の有効性が確認できた．

参考文献
1) 岩澤道徳 : 有限要素法による非線形平衡方程式の数値計算法に関する研究, 長岡技術科学大学大学院
修士論文,2001.2
2) F.Itoh and K.Nogami : On the tracing calculationof equilibrium path for imperfect systems, Proc.
JSCE, Struct.Eng./Earthq.Eng., Vol.3, No.1,1986.4.
3) T. Chaisomphob, W. Kanok-Nukulchai and F.
Nishino : An automatic arc lenght control algorithm for tracing, equilibrium paths of nonlinear structures, Struct. Engng./Earthquake Engng., JSCE, Vol.5, No.1, pp.205-208, 1988.4
4) W. Szyszkowski and J. B. Husband : Curvature
controlled arc-length method, Computational Mechanics, Vol.24, pp.245-257, 1999.