2023-06-14T11:16:14.033山口佳春
日本建築学会構造系論文集 第86巻 第790号, 1607-1614，
2021年12月
J. Struct. Constr. Eng., AIJ, Vol. 86, No.790, 1607-1614, Dec., 2021
DOI https://doi.org/10.3130/aijs.86.1607【カテゴリーⅠ】

不安定性を含む弾塑性解析のための剛性マトリクスの固有値操作
不安定性を含む弾塑性解析のための剛性マトリクスの固有値操作
MANIPULATION OF EIGENVALUES OF A STIFFNESS MATRIXMANIPULATION OF EIGENVALUES OF A STIFFNESS MATRIXFOR AN ELASTO-PLASTIC ANALYSIS INCLUDING INSTABILITY
FOR A ELASTO-PLASTIC ANALYSIS INCLUDING INSTABILITY
昭 夫＊1，江島 ありさ＊2

堀Akio HORI
and Arisa EJIMA
堀 昭夫*1，江島 ありさ*2

Akio HORI and Arisa EJIMA
Eigenvalue manipulation of a stiffness matrix is newly proposed. The method of the eigenvalue manipulation is defined by
Equations 2 and 5 for symmetrical and asymmetrical matrix, respectively. This manipulation magnifies the unbalance forces in
only directions of eigenvectors with negative eigenvalues, and does not affect in the other directions with positive eigenvalues.
In a simple inverted bar model, first order equilibrium leads its results to an unstable balanced path under compression.
Additionally, its stable balanced paths are based on a solution of a cubic equation. Analytical results of 3D cantilever column
model with 4 elastoplastic springs are very successful.

Keywords : Stiffness matrix, Eigenvalue, Manipulation, Elastoplastic, Unstable balanced path
剛性マトリクス，固有値，操作，弾塑性，不安定釣合経路
次以上の剛性)を考える意味が相対的に低下するという立場をとる。

１．はじめに筆者らは過去に，試行的剛性選択法1),2)を提案して，剛性マトリ
3)

そのうえで，不安定釣合経路への落ち込みを回避する目的で，1 次

クス作成時における除荷･載荷の仮定が循環剛性選択過程 に陥った

剛性マトリクス(以下では単に剛性マトリクスと記す)の固有値を操場合でも，試行した仮定の中から不釣合の指標が相対的に小さいも作する解析法を提案する。

のを選択して解析を続行することにより，大規模立体骨組の弾塑性
大変形解析が可能になる事を示した

4),5)。またこの解析法を拡張し２．非常に単純なモデルの不安定釣合経路て提案した立体鋼骨組の高温崩壊解析法により，大規模骨組の高温

本章では後述する固有値操作の提案に先立ち，非常に単純なモデ崩壊解析ができる事を示した6),7) 。しかしこれらの解析も最終時にルを例に，固有値を操作する可能性について調べてみる。具体的には，Fig. 1(a)に示す両端ピンの倒立棒を考える。図中の P1, P2 は外は，釣合解がうまく得られない現象が認められた。
一方，森迫は，座屈荷重を超えて剛性マトリクスの固有値が負と力，d1, d2 は変位，L は材長，EA は断面軸剛性とし，塑性化や慣性なる領域に，整合剛性の下で不安定釣合経路をたどる問題3)があり，力は考えない。倒立棒なので横力 P2 = 0 とし，鉛直変位 d1 を下方

近年の解説書8)においても十分注意しなければならないとしている。

に変位制御する場合を考えて，初期状態で横変位 d2 = 0 とすると，不安定釣合経路の問題は，弾性分岐に関する多くの研究があり，当然ではあるが，倒立棒は軸力 P1 が正なら安定，P1 = 0 で分岐経路

鋼構造座屈設計指針の解説9) や，書籍10) ，総説11) などにも現状がまを持ち，P1 が負なら不安定である。いま，倒立棒の釣合を Fig. 1(b)，2 モード摂の局所座標系で考えてみる。図中に示すように局所座標系方向を添

動展開を行った研究13)，多重固有値を分離した研究14)などもある。

字 x,y で，前ステップ終了時の値を添字 0 で，増分を Δ で示すと，とめられている。また複雑な分岐経路を示した研究12)

また不安定性を含む弾塑性解析のための方法として，動的陽解法(動的緩和法，システムダンピング法，マススケーリング法)の比較
研究15)や一部に質量を入れた飛び移りが起きる架構の解析16) もある
が，質量や減衰の設定により，応力波伝播

15)や動的座屈荷重17)の問題を持ち込んでしまう可能性がありうる。

P 1, d 1 (d 1 control)
P 2, d 2(P 2 = 0)

EA, L

P x, Δd x,
P x 0, d x 0

P y, Δd y,
P y0

L = 1000[mm]θ0

確定しない場合が存在し1),2)，高次の剛性(摂動法で用いるような 2*1 小山工業高等専門学校 教授・博士(工学)＊1
小山工業高等専門学校
教授・博士（工学）*2 ＊
小山工業高等専門学校
専攻科生
2
小山工業高等専門学校 専攻科生

E = 2･10 5[N/mm2]
Δd 1 = -10 -5[mm]

本論では，不安定性を含む弾塑性解析においては，除荷－載荷の剛性仮定の問題があるために 1 次の剛性(接線整合剛性)さえ十分にA = 1000[mm2](a) model(b) local coordinate(c) in analytical cases

Fig. 1 Study model

Prof., Dept. of Architecture, Oyama College, National Institute of Technology, Dr. Eng.
Prof.,of Architecture,
Oyama College,
National
Instituteof Technology,
Grad.Dept.
Student,
Dept. of Architecture,
Oyama
College,
National
Institute ofDr.Eng.
Technology
Grad. Student, Dept. of Architecture, Oyama College, National Institute of Technology— 1607 —2023-06-14T11:16:14.033山口佳春

下式が得られる。

場合の結果を，同図にマゼンタの実線で示す。但し，固有値操作とPx ≒ EA ( ε x0 + Δ ε x + 0.5 Δ ε y2 )(1a)

Py = 0 ≒ Px Δ ε y(1b)

いっても，本章の単純な検討モデルでは，変位制御すると剛性マト
リクスがスカラーになるから，単純に剛性が負であれば -1 を乗じて正値化しただけである。図では 36 ステップ目で d 2 が 0.43 まで成ここに ε x0 , Δ ε x , Δ ε y はそれぞれ，d x 0 , Δd x , Δd y を(初期)材長 L長しているが，これは d 2 初期値の 2 38.6 倍であるから，1-36 のステで除した値とする。なお式(1b)の最初の等号は力の方向を変えない

ップ間において不安定釣合経路からの距離が概ね 2 倍前後に拡大しという境界条件を意味する。式(1a),(1b)より，Δ ε y に関して 3 次，Δ ε xに関して 1 次の，以下の方程式が得られる。

続けた結果として，本来の分岐経路の 1 つへ至ったものと解釈できる。

Δ ε y ( ε x0 + Δ ε x + 0.5 Δ ε y2 ) = 0(1c)

上式は，Δ ε y が 0 の横変位しない解と，括弧内の Δ ε y および Δ ε x に

本章の検討は非常に単純な例ではあるが，剛性マトリクスの負の固有値を正値に操作することで，不安定釣合経路から離してやり，分岐経路の 1 つへと向かわせる事の可能性が示された。

関する 2 次および 1 次の式が 0 になる解があるが，これより以下の３．対称な剛性マトリクスでの固有値操作可能性が示唆される。
①横変位のない基本経路の解と，横変位のある分岐経路の解が独立

前章の解析結果を参考に，対称な剛性マトリクスでの固有値操作的に存在している。このため， 1 次の剛性だけを用いると，分岐として，下式を提案する。なおベクトルを斜体，マトリクスを立体経路を目指さずに，横変位のない安定(Px >0)または不安定(Px <0)で表記する。

な基本経路上の点を目指してしまう

注1)。

②このような基本経路へと向かう問題は，森迫らが詳述した，座屈
荷重を超えた領域において整合剛性の下で不安定釣合経路をたど

mK = K - Σ mi λ i u i u i
i

(2)

T

ここに，mK は操作後の剛性マトリクス，K は操作前の(通常の)剛性マトリクス，λ i および u i は K の i 番目の固有値および固有ベク

3)る問題 と部分的に似た状況にあると思われる。
③これらをもとに，剛性マトリクスの固有値が負となった場合には，

トル(ノルムは 1 に基準化)，Σi は負の固有値と対応する部分の総和，

1 次の剛性に従うと不安定な経路を目指してしまうと仮定する。

m i は後述する固有値操作の定数，とする。上式による mK は以下の

この場合，1 次の釣合式による解は，力の釣合は満足するものの特性を持つ。

不安定であるから，その解を指向しても実際の挙動からは遠ざか①対称な K に対して，線形代数の定理より固有ベクトルの直交性が成立するから，式(2)の両辺に後から j 次の固有ベクトル u j を乗ることになる。
④このため本論では，不安定釣合経路に向かわないよう，剛性マト
リクスの負の固有値を正値へと操作する事を考え，それによる挙動を調べてゆく。
Fig. 1 の倒立棒の数値解析結果を Fig. 2 に示す。なお解析用の条
件は Fig. 1(c)によっており，変位制御増分の Δd1 を-10 mm とし-5たが，d 1 の初期値も-10 -5mm にして初期軸力をわずかに圧縮側とじれば，mK u j = K u j， 但し i ≠ j の場合(3a)となる。つまり，固有値操作をしない他の固有ベクトルに関わる
変位増分については，K を mK に代えても影響が起きない。
②負固有値を持つ i 番目の固有ベクトルに関わる変位増分に対してした。Fig. 2 において，◆付の青線および○付の赤線は d 2 の初期値

は，剛性マトリクスの固有値が(1 - m i)倍される。これは，式(2)がそれぞれ 0.04mm および 0.05mm の場合を示すが，いずれも固有の両辺に後から u i を乗じて考えると(但し負固有値が 1 つだけの値操作は行わずに，鉛直下方に d 1 を変位制御した計 50 ステップ程

場合を略記するが)，度の解析の結果である。なお各ステップでの収束計算は行わず，次のステップで不釣合解消力を作用させている注2)。同図を見ると，2mK u i = K u i - mi λ i u i= (1 - m i) λ i u i(3b)つの場合とも初期の数ステップの振動の後に，圧縮軸力下の不安定となるから，容易に確かめることができる。この定数 m i の値に釣合経路をたどる場合(d 2 = 0.04)と，軸力 0 の分岐経路をたどる場

ついては，相当に複雑な剛性選択を含む場合の結果も勘案して定合(d 2 = 0.05)に分かれた様子がわかる。

めるべきであるが，本論の範囲では特に言及のない限り m i を仮

一方，d 2 にわずかな初期値(10-12)を与えて，固有値操作を行った

値の 2.0 として検討を進める。この場合，K の負固有値を-1.0 倍したものが，mK の正固有値へと操作される。
③上記①②より当然ではあるが，mK の固有ベクトルは K のそれとinitial d 2 = 0.05

変わらない。また対称な K に対して，mK も対称となるが一般には密行列になる。

initial d2 -->0
manipulated④負の固有値が多数ある場合でも，上記の固有値操作には特段の支障が起きないと考えられる。この点は，弾性分岐問題において 0
固有値が複数ある場合に工夫する方法13)や，弾塑性分岐問題におinitial d 2 = 0.04

3)

いて負固有値の固有ベクトルの正負両方向に試行する方法 ，にはない特長となる可能性が考えられる。但し多重固有値に対してFig. 2 Results of inverted bar— 1608 —は mi の値を同一にすべきと想像される。

2023-06-14T11:16:14.033山口佳春

Force

along
(1-mi)λi

unstablebalance point

Force

x, Δx

unstablebalance point

Disp.

L = 1000[mm]e y = ±20[mm]e z = ±20[mm]k e = 300[N/mm]
k p = 3[N/mm]
P y = 1.5[N]d y = 5･10-3[mm]F y = F z = 4.8･10-4[N]Δx = -4･10-4[mm]y, F yDisp.

unbalance force
magnified by mi*z, F zLPy(e y, e z)magnifiedby mi*kpkedyalong λi(a) linear(b) quadratic
Fig. 3 Magnification of unbalance force(b) bi-linear spring(a) model(c) in analysis
Fig. 4 Spring column model

なお当該ステップの前後における剛性変化が小さい場合には，式(2)を用いて負固有値を(1 - m i)倍する事による当該固有ベクトル方向の不釣合力の拡大縮小率 m i* は，ステップ内での収束計算を特にしなければ，
m i* = 1.0 -1.0/(1 - m i) = m i /(m i - 1 ) ,if linearA

Bとなるはずである(Fig. 3(a))。これより m i が 2.0 の場合には m i*
も 2.0 になるが，もしも不釣合力をより積極的に拡大したい場合に

B

A

(4a)
(a) axial force

は上式から m i を設定する事もできよう。一方，もし剛性変化が大

pp
ep

ee
ee

きく，荷重―変位関係が Fig. 3(b)のような放物線の場合には，

pp
pp

m i* = {1.0 -0.5/(1 - m i)}2 = (m i - 0.5 )2/(m i - 1 )2 , if quadratic(4b)(b) lateral disp.

Aとなるものの，m i が 2.0 の場合の m i* は 2.25 にしかならず，Fig. 3(c) y,z disp.

B
pp
ee

(λ2 negative)(d) eigenvalues(a)の場合の 2.0 との相違は大きなものではない。このため，m i の
仮値 2.0 は，滑らかな剛性変化には影響されにくい範囲にあるよう+zに思われる。
+y４．単純な弾塑性柱モデルでの解析結果
前章の固有値操作を Fig. 4 のモデルに対する解析に適用して，そAB(e) angle of eigenvectorの効果を調べる。同モデルは，剛体棒の脚部に 4 つのバイリニアばねを持つ片持柱立体モデルで，圧縮降伏軸力(6N)に達した瞬間に，
塑性座屈軸力(強軸弱軸がなく 2 方向とも 4.8N)を超過するモデルとなっている。もし固有値操作をしない場合には，整合剛性下で不安
定釣合経路をたどる可能性がある

3),8)

モデルである。なおモデルは特定の柱材を想定したものではなく，形状変化をあまり大きくしな(f) each spring
Fig. 5 Results of spring column model

いように降伏変位を小さめに設定している。
このモデルに対し，軸方向(-x 方向)に変位制御して 100 ステップの解析を行った。なお変位制御増分は圧縮降伏変位の 1/12.5 とした。

向(+y,+z 方向)であったものが，途中の B 点で 90° 方向(+z 方向)に

また頂部の水平 2 方向(+y,+z 方向)には，それぞれ 4.8･10-4N の荷変化している。剛性マトリクス固有値の操作前の値を fig.(d)に示す。

重不整を与えた。この不整は，塑性剛性下では y,z 方向とも頂部変

但し fig.(d)の菱形内の e および p は，各ばねが弾性(除荷)剛性か塑

位 0.1mm の増加に対応し，弾性剛性下では y,z 方向ともばね降伏力性剛性かを示し，Fig. 4(a)の 4 つのばね位置と対応するように位置の± 0.4%だけばね力を増減させる曲げに対応する。なお本論の解析をずらして描いてある。推移を以下にまとめる。

例では慣用の試行修正過程3) を用いており，具体的には，固有値操i)初期には固有値が 2 つとも正であるが，ii) 4 つのばねが降伏す

作後に求めた各ばねの変形増分が，剛性マトリクス作成時の除荷･載ると 2 つの固有値とも負になって固有値操作の対象となり，その結荷の仮定と矛盾する場合には，矛盾したばねの除荷･載荷の仮定を変果，循環剛性選択過程を経ずに，慣用の試行修正過程 3)により図中更して剛性マトリクスを作り直し，変位増分を再計算している。
解析結果を Fig. 5 に示す。まず圧縮力―圧縮変位関係(fig.(a))は図中 A 点で降伏するバイリニアとなっている。次に，ばね降伏後(A

左下(-y,-z 方向)のばねが除荷された(A 点)。なお操作前の 2 つの負固有値は，少なくとも有効数字 7 桁が一致して多重固有値に近かった。

点以降)は横変位の増加が見られ(fig.(b))，塑性座屈した事がわかる。

iii)その降伏状態がしばらく続いたが(A-B 間)，その間の固有値

ところが fig.(c)の横変位方向(座屈方向)を見ると，初期には 45° 方

は 1 つだけ負であり，固有値操作によってその固有ベクトル方向の— 1609 —2023-06-14T11:16:14.033山口佳春

不釣合力が拡大を続けて，44 ステップ後に右下のばね(+y,-z 方向)

負値(-y,-z 方向)に限定されていて，固有ベクトル方向も 45° 方向かも除荷された(B 点)。但し 2を考えると，倍精度の数値

ら変化しなかった。fig.(a)に 45° 方向変位を示すが，ばね降伏後(A

計算誤差が発端となって初期不整のない方向に座屈方向が変化した点以降)の横変位は荷重不整と逆向きで，かつ C 点以降は座屈変位と推察できよう。iv) 2 つのばねが除荷された後は，2 つの固有値とが減少をしている。fig.(b)に圧縮力－ 45° 方向変位の関係を示すが，も正になっている。

C 点以降では圧縮力が増加しつつ座屈変位が減少するという，森迫

44

≒ 10

13.2

固有ベクトルの方向を fig.(e)に示す。固有ベクトルの方向は 45°
方向(荷重不整と同方向)が続いたが，2 つめの除荷発生(B 点)以後

ら 3) が詳述した不安定釣合経路に陥ったと判断される。fig.( c)に剛性マトリクスの固有値を示す。推移は以下のようであった。

は方向が変化している。なお A 点では負固有値が多重に近いが，fig.

i) 4 つのばねが降伏した直後に循環剛性選択過程に陥った。1 つ(e)を見る限り，A 点の前後で固有ベクトル方向は変化していない。

の剛性仮定は，4 つのばねとも塑性継続(固有値が 2 つとも負)，もfig.(f)に各ばねの力と変位の関係を重ね書きして示す。なお Fig. 4

う 1 つの剛性仮定は，右上のばねだけが除荷発生(固有値の 1 つが( a)の各ばね位置に合わせて，それぞれ初期位置をずらして描いて負)，となっていた(A 点)。このため本例では，後者の剛性仮定を採ある。座屈発生後も左下のばねはほとんど塑性変形していない(初期用して解析を続行した(なお前者の仮定を採用すると解析結果が発散

位置が原点の曲線)。

的になった。また慣用の試行修正過程3) を用いたため，これら以外

なお本解析対象の性質(不整がない場合)は，以下のように想像さの剛性仮定は試行していない)。

れる。4 つのばねが降伏すると，鉛直方向に変位制御した関係で，ii)この循環剛性選択過程が C 点まで続いたが，C 点以後は 4 つの水平面内の 0° ～ 360° の全方向に不安定になるが，除荷･載荷に関しばねが塑性継続で 2 つの固有値が負のままで，剛性仮定と変形増分て停留(除荷)のばねが 1 つ発生した状態では不安定性が残り，停留の方向が整合する釣合経路をたどっている。しかしこれは，固有値(除荷)のばねが 2 つ発生して安定化するように見受けられる。本章の符号から考えて，2 方向に不安定な釣合経路である。

の解析結果では，4 つのばねの降伏直後に 45° 方向の初期不整によ

fig.(d)に各ばねの力と変位の関係を重ね書きして示す。4 つのばって 1 つのばねだけが停留的な挙動になったが，固有値操作で不釣ねが塑性を継続しているが， fig. 5( f)とは異なり，左下のばねの塑合力が拡大した結果，2 つのばねが停留的になる方向へと移行した。

性変形が最も大きい。

なお 2 つのばねが停留的になる変形方向(fig.(c)の B 点以降)は固有

以上のように，固有値操作をしない本例の解析結果では，①本来ベクトルの方向(fig.(e)の B 点以降)とは一致していない。また解析とは逆方向に座屈した後に，② 2 方向に不安定な釣合経路をたどる，

結果の B 点以降の変形は 90° 方向であったが，90° 方向か 0° 方向かという結果に陥った。

は本解析例では偶然性(数値解析誤差ほか)によると思われる。
６．降伏力を10倍または75倍にした場合の解析結果
Fig. 5 の解析に対し，ばねの降伏力･降伏変位および変位制御増分

５．固有値操作をしない場合の解析結果
Fig. 5 の解析に対し，もしも固有値操作をしない場合の結果をFig. 6 に示す。最初に，図は省略するが，圧縮力―圧縮変位関係はFig. 5(a)と全く同様であった。但し，座屈後の横変位は 45° 方向のΔx をいずれも 10 倍にした解析を行った。結果を Fig. 7 に示し，元の結果(Fig. 5)との相違点について述べる。
まず圧縮力―圧縮変位関係(fig.(a))では，座屈後の耐力低下が見られる。横変位は 4 つのばねが降伏した A 点以降で増加している

A

A

C

(fig.(b))。しかし A 点の解析ステップでは固有値操作により直ちに除荷発生するのではなく，A-D の 4 ステップの間は 2 つの固有値が負のままで，D 点で初めて左下のばねが除荷されて(fig.(d))，負固C

有値の数が 1 つに減じた。但し A-D 間では，1 ステップ毎に横変位増分が倍増する現象が起きていた。A 点から 32 ステップ後の E 点(a) +y+z disp.
eeee(b) axi. force vs +y+z disp.

pepp (λ1 positive?, negative? )

A
pp
pp

( λ2 negative)Cで(右下ではなく)左上のばねが除荷されて(fig.(d))，初期不整のない方向に座屈方向が変化した(fig.(c))。
次に，ばねの降伏力･降伏変位および変位制御増分 Δx をいずれもpp
pp(λ1, λ2 negative)(c) eigenvalues

75 倍にした解析を行った。なお，もし 75 倍でなく 80 倍にすると，不整を無視した弾性座屈荷重( 480N)が圧縮降伏荷重と等しくなって，限界細長比と対応してくる。また本例の解析では横変位増分が過大になるステップが見られたので，横変位増分を 10mm 以下に制限している。この点は，本論の固有値操作が 0 固有値に対しては何ら影響を持ち得ない事から，解析対象によっては変位増分の上限を設ける等の必要性を意味するであろう。
解析結果を Fig. 8 に示し，前例の結果(Fig. 7 )との相違点について述べる。但し横変位方向および 1 つの固有ベクトル方向は 45° 方(d) each spring
Fig. 6 Results without eigenvalue manipulation— 1610 —向から変化しなかったので，それらの結果を省略する。まず圧縮力―圧縮変位関係(fig.(a))では，座屈後に大きな耐力低下(F 点-G 点)

2023-06-14T11:16:14.033山口佳春が見られるが，H 点-I 点でも耐力低下が見られる。座屈後の G 点では片持柱の傾きが 2%近い(fig.(b))。固有値は，A-F 間では 2 つとADも負であったが，G 点以降は 2 つとも正になっていて，G-H 間では 3E

つのばねが除荷状態となっている(fig.(c))。また H 点の直後で，左E(a) axial force

A

下のばねが引張降伏し，かつ右下および左上の 2 つのばねが再降伏(b) lateral disp.

つのばねが除荷されて，固有値は 2 つとも正になった(fig.(c))。各(圧縮)して，2 つの固有値が一度負になったものの，次の I 点で 2ばねの力と変位の関係を見ると，左下のばねが引張側(負側)に大きく塑性変形している(fig.(d))。

E

解析対象としたモデルで，不整を無視し，かつ，ばねの力と変位の関係が弾性勾配上，圧縮耐力直線上もしくは引張耐力直線上にあ(c) y,z disp.

ると仮定して求めた(複数の)釣合曲線を fig.(e)に示す。マゼンタ線，(f) each spring

ee
ee

ep
ep

pp
ep

D

ep

pp

ep

pp

それぞれ 1･3･0(e e )，3･1･0(e p )，1･2･1(pe )，3･0･1(pp )を仮定している。また各線とも，構成則と整合しない区間は点線で描いた。
i) 同図を見ると，A 点からの釣合曲線は，3 つのばねが除荷され
ep

E(λ2 negative )A

青線，緑線，黒線は圧縮降伏･弾性･引張降伏のばね数の仮定として，るマゼンタ線(e e )であり，軸変位の戻りを伴っている。図中には解析結果(赤線)も重ね書きしたが，解析結果では A 点から釣合解のな

pp
pp

い方向に変位制御をしていて，不釣合力が大きくなった 6-7 ステッ(d) eigenvalues

pp

ep

プ後の座屈により大きく耐力低下して，( e p )を経て(e e )に変化し+zてマゼンタ線に沿っている。
ii) しかしその直後に，マゼンタ線と青線が交差していて，釣合

E

pp

ep

曲線が青線(e p )に移っているが，解析結果は(e e )のまま若干推移

+y

している。
(e) angle of eigenvector
Fig. 7 Results using strengthened springs to 10 timesiii) ところが図中の h 点付近では，マゼンタ･青･緑･黒の 4 つのppep

線が密集していて，かつ釣合曲線が青線( e p )から緑線( pe )に移っていて，3 つのばねの状態変化も示唆される。解析結果では 8 ステ
pp

ep

ップ程度遅れて(pp )を経て緑線の降伏状態(p e )になって耐力低下があり，以後は緑線へと漸近している。

A

F
G

以上に示した Figs. 5,7,8 の片持柱モデルの 3 つの解析例では，固
H
I

G

I

有値操作を行った結果，循環剛性選択過程を経る事なく，初期にはH

荷重不整のある 45° 方向へ座屈し(Figs. 5,7,8)，数 10 ステップの後には初期不整のない方向へ座屈方向が変化したり(Figs. 5,7)，釣合(a) axial force

G
A

pp
pp

F

pp
ep

(b) lateral disp.

ep
ee H

解のない方向への変位制御であっても数ステップ遅れるものの座屈をしたり( Fig. 8)，これらの経過の後にいずれも負の固有値が解消されてゆく結果が示された。

ep
pe

I

７．非対称な剛性マトリクスでの固有値操作

pp
pp

以上では剛性マトリクスの対称性を暗黙的に利用してきたが，立体で回転自由度があると，座標更新の方法にもよるが幾何剛性マト(c) eigenvalues

リクスが非対称になる場合がある注3)。非対称マトリクスは，固有値が複素数になる可能性があり，固有ベクトルの直交性が一般に成立Aしなくなる注4)。本論では，非対称な場合も含む剛性マトリクスの固hH(d) each spring(e) ideal curves
Fig. 8 Results using strengthened springs to 75 times

有値操作として，下式を提案する。
mK as = K - U j i m i λ i U ji

-1

(5)

ここに，mK as は操作後の(対称または非対称な)剛性マトリクス，Kは操作前の(対称または非対称な)剛性マトリクス，λ i は K の固有値を並べた対角マトリクス，m i は m i を並べた対角マトリクス，U j i
は K の固有マトリクス(ノルムを 1 に基準化した i 番目の固有ベク— 1611 —2023-06-14T11:16:14.033山口佳春

トル u i を i 列に配した行列)，とする。なお添字 i は i 番目の固有値領域を通過した後にもいくらかの荷重低下を示した。このように本･固有ベクトルに対応する事を示す。式(5)の mK as は以下の特性を論の弾塑性解析結果からは，固有値操作が不安定性を有する問題に持つ。

対して，数ステップ(以上)の遅れを伴うが非常に有効である事が例①いま系の増分変位ベクトル Δx を，各固有ベクトル方向に分解し証された。

た場合の係数ベクトル α i を用いて，次式で表す。

本論では，固有値操作の定数 m i のより適切な値については述べΔx = U ji α i(6)

このとき，式(5)の両辺に後から Δx を乗じると，
m K as Δx = K Δx - U j i m i λ iていない。また実骨組レベルの解析結果や剛性マトリクスが非対称な場合の解析結果も示していない。これらについては今後の課題である。また大規模な骨組への適用を考えると，大規模行列の固有値･

α i = K Δx - U j i {m i λ i α i}(7)

固有ベクトルの計算コストが高いため，サブストラクチャ法を適用して，局部的に固有値操作する方法が考えられる。この観点からは，となる。ここに {m i λ i α i} は，スカラー m i , λ i , α i の 3 つによる本論の固有値操作が，あくまでも連立 1 次方程式の求解前の剛性を積を i 番目の成分とするベクトルとする。上式を式(3b)と比べて操作している点も重要と思われる。仮に連立 1 次方程式の求解後の考えると，mK as は，i 番目の固有ベクトルに関わる変位増分に対変位増分から負固有値に対応する固有ベクトル方向の変位増分を拡して，剛性マトリクスの固有値が(1 - m i)倍されている事がわかる。

大縮小する手続きであったとしたら，サブストラクチャ法と併用しなお，mK as の固有ベクトルは K のそれと変わらない。

づらくなるからである注5)。なお本論の固有値操作を弾性問題に適用②以上により，非対称な剛性マトリクスに対しても，式( 5)により固有値操作が支障なく定式化できる事がわかる。また，もし K が対称である場合には，線形代数の定理より U ji

-1

は U ji と一致す
注

③式(5)の固有値操作は，固有ベクトルの直交性を利用していない。
また(非対称性が弱く用価値があるかもしれない。

Tるから，式(5)は式(2)を包含した表現である事もわかる。
注4)

する事は排除しない。重解に近い負固有値を持ちうる問題などで適)固有値が実数であれば，固有値操作の物
注4)

理的な意味も明瞭である。しかしながら，(非対称性が強く)固

有値が複素数になった場合には， I)線形代数の定理より複素共役となる固有値の正負の判定をどうするか，II)線形代数の定理より
複素数となる固有ベクトル方向の変位や力を物理的にどう考えるか，等の問題が考えられよう。
８．結語
本論では，剛性マトリクスの固有値を操作する方法を新たに提案した。これを固有値操作と呼ぶ。なおこのような操作が理工学の他

注 1) 式(1a),(1b)から局所座標系での剛性マトリクスを記述すると，
ΔP x
ΔP y

≒ EA

0.5 Δ ε yε x0 + 0.5 Δ ε y2

1
Δεy

Δεx

(n1.1)

Δεyと書けるかもしれない。ここに，上式中の下線部は増分変位に関し 2 次の剛性の項，下線部は増分変位に関し 3 次の剛性の項，になっている。後述
する本論の解析では，1 次の剛性(無下線部に対応)のみを使って増分計算しており，2 章の解析では上式の無下線部を全体座標系に変換した後に，d1を変位制御する条件下での剛性(この場合はスカラーになる)の固有値を場合により操作して，仮に操作しない場合には P 2 = 0 を目指すような解を求めている。
なお上式の無下線部のみを用い，かつ固有値操作しない場合には，
ΔP 1
ΔP 2

≒

EA 1 -θ0
L θ0 1

1
0

0

≒

EA 1
L θ0

Δd1
Δd2

ε x0

1 θ0
-θ0 1

Δd 1
Δd 2

分野で行われているか否かは，著者らによるインターネット検索ではわからなかった。固有値操作の式は，対称マトリクスに対しては
式(2)，非対称マトリクスに対しては式(5)，で与えられる。

θ0

ε x0

(n1.2)

固有値操作の提案に至った経緯は，Fig. 1(a)のような非常に単純とできるから，これに境界条件の P2 = 0 を Δd 1 → 0 の条件下で適用すると，な系において，釣合解が 3 次方程式の解になり，もし 1 次方程式の(n1.3)

EA･ε x0･θ0 + ΔP 2 = 0 = EA･ε x0 (θ0 + Δd 2 /L)

解を用いると不安定な釣合経路解を指向してしまう，という点にあとなって，横変位 d 2 を 0 にする Δd 2 が解になるから，軸力の正負に関係る。このような場合の 1 次方程式は，工学的利用の観点からは力をなく基本経路上の点を目指す現象が起きてしまう。
注 2) 本論の解析では，各ステップでの変位制御増分を一定にし，ステップ内での収束計算を行わずに，当該ステップで生じている不釣合力を次ステッ
プの外力に逆符号で加えている。本論では固有値操作を新たに提案する観点から，解析手法をなるべく単純にしたままで，固有値操作による解析結果がどうなるかを調べてゆく。このため，増分解析で潜在する他の問題，

釣り合わせる意味がそもそも薄い事になる。本論の固有値操作は，
このような不安定釣合経路からは遠ざかるように，負固有値に対応する固有ベクトル方向の不釣合力を拡大してゆく目的･性質を持つ。
不安定釣合経路から遠ざける操作が，直ちに安定釣合経路に落ち着く事を意味するわけではない。しかしながら，Figs. 5-8 の弾塑性解析結果を通じて示された事は，①多重に近い 2 つの負固有値があったが循環剛性選択過程を経る事なく座屈して，かつ座屈方向が途中から不整のない方向へと変化した(Fig. 5)。②固有値操作をしない
解析結果では循環剛性選択過程に陥り，除荷剛性を仮定して解析を続けると，本来とは逆方向に座屈した後に 2 方向に不安定な釣合経路をたどった( Fig. 6)。③限界細長比に近いような系で，釣合解のない方向に変位制御した解析では，不釣合力が成長した後に大きな荷重低下を伴って座屈して，不整のない場合の釣合曲線が密集する— 1612 —例えば，不安定性がある場合に解析中の制御増分の大きさを変更すべきか否か，またもし変更する場合には増分変位や不釣合力の大きさから一義的に決めてよいか，等の問題には立ち入らない。これらの点については，より複雑な弾塑性問題で除荷･載荷の整合剛性が容易には定め難い場合 例えば ,1),2)，などへの適用性も含めて検討する必要があるかもしれない。またス
テップ内での収束計算については，著者らは，除荷･載荷の接線剛性が確定しない(矛盾する)場合には収束しにくい場合がありうると考えている。
注 3) 右手系で例えば，曲げモーメント外力 M = (Mx, M y, Mz)が Δθ = (Δθx,Δθy, Δθz) だけ回転する場合の幾何剛性 MKG は，

MKG =

- Mz

0

Mz
0

- My

My

-Mx

0

Mx

,

但し，ΔM = MKG Δθ ，とする

(n3.1)

2023-06-14T11:16:14.033山口佳春となって，交代行列になる。なお線形代数の定理より，実交代行列の固有値は全て純虚数または 0 である。
注 4) 下記の非常に単純な実マトリクスB=1q

森迫清貴ほか: 建築構造における強非線形問題への数値解析による挑戦, 2.1
節, 応用力学シリーズ 13, AIJ, pp.7-15, 2018p1

(n4.1)の固有値 λ および固有ベクトル u 1, u 2 は，
λ = 1 ± (p q )[ u 1, u 2 ] =(n4.2)

0.5pp(p q)0.5- (p q)0.5

Nonlinear Problems in Building Structures, Section 2.1, Applied

Mechanics Series 13, AIJ, pp.7-15, 2018 (in Japanese),

但し p q ≠ 0(n4.3)とできる。つまり B では，固有値･固有ベクトルが実数になる条件は，p q≧ 0，つまり非対角要素の符号が異ならなければよいので，①非対称性の程度が大きくなければ固有値が実数になる。
しかし一方，u 1 と u 2 の交角の余弦は，(p-q)/(p+q) となるので，p ≠ q
つまり，②非対称マトリクスであれば固有ベクトルが直交しなくなる。
よって①②より，B では固有ベクトルの直交条件の方が固有値･固有ベク
トルの実数条件よりも厳しくなっている。ただし B は p q = 0 の場合には，多重固有値になって固有ベクトルの方向に任意性が生じる(p,q とも 0 の場
合)かまたは，固有ベクトルが線形独立でなくなる(p,q の片方だけが 0 の場合)。
注 5) 但し逆に，固有値操作後の剛性が密行列になる問題を回避する観点から，固有値操作せずに変位増分を求めて，負固有値に対応する成分を拡大縮小する，という方法も想定可能である。しかし固有値･固有ベクトルの計算コストはこの場合も同じである。

参考文献
1) Hori, A. and Sasagawa, A.: A Combined Non-linear Analytical Method

for a Space Frame using one Dimensional Finite Element Method,
Journal of Structural and Construction Engineering (Transactions of

AIJ), No. 490, pp.139-147, 1996.12 (in Japanese)

堀昭夫, 笹川明: 試行的な剛性選択法を用いた 1 次元有限要素法による立体
骨組の複合非線形解析法, 日本建築学会構造系論文集, 第490号, pp.139-147,
1996.122) Hori, A. and Sasagawa, A.: Large Deformation of Inelastic Large Space
Frame. I: Analysis, Journal of Structural Engineering, ASCE, Vol. 126,
No. 5, pp.580-588, 2000.5
3) Uetani, K., Nakamura T., Morisako, K. and Ishida, S.: A Method ofGenerating Consistent Stiffness Matrices for Incremental Analysis of
Critical Behavior of Eelastic-plastic Structures, Journal of Structural

and Construction Engineering (Transactions of AIJ), No. 445, pp.67-78,
1993.3 (in Japanese)
上谷宏二, 中村恒善, 森迫清貴, 石田修三: 弾塑性構造物の臨界挙動解析のための 整合剛性行列形成法 , 日本建築学会構造系論文報告集 , 第 445 号 ,pp.67-78, 1993.3
4) Hori, A. and Sasagawa, A.: Combined Non-linear Analyses of a Large

9) AIJ Recommendations for Stability Design of Steel Structures, 4th
Edition, Chapter 10, 2018 (in Japanese)
日本建築学会: 鋼構造座屈設計指針, 第 4 版, 10 章, 201810) Neto, EA., Perić, D. and Owen, DRJ.: Computational Methods for
Plasticity (trans supervisory) Terada, K., Morikita Publishing, 2012
Neto, EA., Perić, D. and Owen, DRJ.: 非 線形有限要素法, 寺田賢二郎監訳,
森北出版, 201211) Champneys, AR. et al.: Happy Catastrophe: Recent Progress inAnalysis and Exploitation of Elastic Instability, Frontiers in Applied

Mathematics and Statistics, Vol. 5, Article 34, doi: 10.3389/fams.2019.
00034, pp.1-30, 2019.7

12) Obiya, H. et al.: A Study on Equilibrium Shape of Tensegrity

Structures with Virtual Stiffness through Some Numerical
Experiments, Journal of JSCE, Ser. A2 Applied Mechanics, JSCE, Vol.
68, No. 2, pp.45-56, 2012 (in Japanese)

帯屋洋之ほか4名: テンセグリティ仮想剛性構造の釣合解特性に関する数値
実験的研究, 土木学会論文集A2（応用力学）, 第68巻, 第2号, pp.45-56, 2012

13) Fujii, F. et al.: Validating a 2-mode Asymptotic Expansion in
Computational Bifurcation Theory, Transactions of the JSME, JSME,
Vol. 81, No. 830, pp.15-00419, 2015 (in Japanese)

藤井文夫ほか3名: 計算分岐理論における2-mode漸近展開の検証, 日本機械学会論文集, 第81巻, 第830号, pp.15-00419, 201514) Saiki, I. et al.: An Improvement of the Stiffness Modification

Method to Find Bifurcated Paths at a Multiple Bifurcation
Point, Journal of JSCE, Ser. A, JSCE, Vol. 62, No. 4, pp.782-793,

2006.10 (in Japanese)
斉木功ほか3名: 多重分岐点における分岐経路探査のための修正剛性法の改良, 土木学会論文集A, 第62巻, 第4号, pp.782-793, 2006.10
15) Tsuda, M. and Hagiwara, I.: Dynamic-Explicit Finite Element Analysis
Methods for Large-deformation Quasi-Static Problems, Transactions of
the JSME, Ser. A, JSME, Vol. 64, No. 622, pp.114-121, 1998.6 (inJapanese)

津田政明, 萩原一郎: 準静的大変形問題の動的陽解法有限要素法に関する基
礎検討, 日本機械学会論文集A, 第64巻, 第622号, pp.114-121, 1998.6

16) Ssuzuki, H. et al.: Stabilities of Steel Frames Subjected to Fire,
Journal of Structural and Construction Engineering (Transactions of
AIJ), No. 571, pp.161-168, 2003.9 (in Japanese)

鈴木弘之ほか2名: 火災加熱を受ける鋼架構の構造安定性, 日本建築学会構造系論文集, 第571号, pp.161-168, 2003.9
17) Timoshenko, SP. and Gere, JM.: Theory of Elastic Stability (trans)
Hasegawa, M., Brain Books Publishing, 1974
Timoshenko, SP. and Gere, JM.: 弾 性安定の理論, 長谷川節訳, ブレイン図書出版, 1974

Space Frame using One-dimensional Finite Element Method, Journal of
Structural and Construction Engineering (Transactions of AIJ), No. 494,
pp. 75-81, 1997.4 (in Japanese)
堀昭夫, 笹川明: 1 次元有限要素法による大規模立体骨組の複合非線形解析,

日本建築学会構造系論文集, 第494号, pp.75-81, 1997.4
5) Hori, A. and Sasagawa, A.: Large Deformation of Inelastic Large Space
Frame. II: Application, Journal of Structural Engineering, ASCE, Vol.
126, No. 5, pp.589-595, 2000.5
6) Hori, A.: Analytical Method of High Temperature Collapse of 3D SteelFrame, Journal of Environmental Engineering (Transactions of AIJ),
No. 586, pp. 9-16, 2004.12 (in Japanese)
堀昭夫: 立体鋼骨組の高温崩壊解析法, 日本建築学会環境系論文集, 第586号,
pp.9-16, 2004.12
7) Hori, A.: Analytical Method for High Temperature Collapse of a 3D
Steel Frame, Fire Science and Technology, Center for Fire Science andTechnology, Tokyo University of Science, Vol. 23, No. 3, pp.208-221,
2004
8) Morisako, K. et al.: Challenge with Numerical Analyses for Strong— 1613 —2023-06-14T11:16:14.033山口佳春

MANIPULATION OF EIGENVALUES OF A STIFFNESS MATRIXMANIPULATION OF EIGENVALUES OF A STIFFNESS MATRIXFOR AN ELASTO-PLASTIC ANALYSIS INCLUDING INSTABILITY
FOR A ELASTO-PLASTIC ANALYSIS INCLUDING INSTABILITY
Akio HORI ＊1 and Arisa EJIMA ＊2

Akio HORI*1 and Arisa EJIMA*2＊1

Prof., Dept. of Architecture, Oyama College, National Institute of Technology, Dr.Eng.
Grad. Student, Dept. of Architecture, Oyama College, National Institute of Technology*2 Grad. Student, Dept. of Architecture, Oyama College, National Institute of Technology2*1＊Prof.,Dept. of Architecture, Oyama College, National Institute of Technology, Dr. Eng.

Eigenvalue manipulation of a stiffness matrix is newly proposed. The method of the eigenvalue manipulation isdefined by Equations 2 and 5 for symmetrical and asymmetrical matrix, respectively.
m K = K - Σ mi λ i u i u i
i

(2)

T

m K as = K - U j i m i λ i U j i

-1

(5)

Herein, for Eq. 2, mK and K are stiffness matrices after and before manipulation, respectively. λ i and u i are i-th
eigenvalue of K and its eigenvector normalized to 1.0. Symbol Σ with i means partial summing within concerning
negative eigenvalues. m i is manipulation constant and set as 2.0 tentatively in this paper. In Eq. 5, U ji iseigenmatrix of K whose i-th column is u i. m i and λ i are diagonal matrices as a list of m i or λ i, respectively. The
eigenvalue manipulation magnifies the unbalance forces in only directions of eigenvectors with negative eigenvalues,
and does not affect in the other directions with positive eigenvalues.
In elastoplastic analyses including instability, stiffness matrix sometimes has difficulty to be determined becauseof how to make a consistent assumption of loading or unloading for every fiber, section or member stiffness.
Moreover, these analyses sometimes fail into an unstable balanced path, which is shown even in a simple inverted
bar model in Fig. 1. In this model, first order equilibrium leads its results to an unstable balanced path in caseunder compression. Additionally, its stable balanced paths are based on a solution of a cubic equation. Therefore, we
proposed the manipulation of stiffness negative eigenvalues to positive ones to avoid falling into unstable balanced
paths. Analytical results of the above simple bar model is successful as shown in Fig. 2, although the manipulation
is not proved toward a stable balanced path for all cases.
More detailed cases are studied using 3D cantilever column model as with 4 elastoplastic springs on the base. The
results using the eigenvalue manipulation are shown in Figs. 5, 7 and 8 and did not fall into a cyclic stiffness
selection process. The results show, i) buckling initially toward 45 degrees parallel to load imperfection (Figs. 5, 7
and 8), ii) shifting buckling direction toward one with no imperfection after several tens of steps (Figs. 5 and 7), iii)
buckling with delay of several steps against displacement control direction of no equilibrium solution (Fig. 8), iv)
two negative eigenvalues in the early stage turning to positive ones after above those processes.
Additionally, results with no eigenvalue manipulation (Fig. 6) fell into a cyclic stiffness selection process andbuckled to the opposite direction against its initial imperfection if selecting the unloaded stiffness assumption. The
results finally failed into a two-directional unstable balanced path with a decrease of buckled deformations.
The eigenvalue manipulation proposed here seems to have broad applicability against instability.
（2021 年 3 月 26 日原稿受理，2021 年 8 月 12 日採用決定）— 1614 —