第一工業大学研究報告第32号（2020）pp.11-16

2 質点系せん断モデルに弾塑性復元力特性を組み込んだ系で観測されるカオス的性質
山田 猛矢 1 ，福永 知哉 2
1

第一工業大学2

工学部

第一工業大学

情報電子システム工学科（〒 899-4395 鹿児島県霧島市国分中央 1-10-2）
E-mail: t-yamada@daiichi-koudai.ac.jp
共通教育センター（〒 899-4395 鹿児島県霧島市国分中央 1-10-2）
E-mail: t-fukunaga@daiichi-koudai.ac.jp

Chaotic Properties Observed in 2-Mass Shear System
Incorporating Elasto-Plastic Force-Displacement Relation
Takeshi YAMADA1 ，Tomoya FUKUNAGA2
1

Department of Informatics and Electronics, Daiichi Institute of Technology
2
Common Education Center, Daiichi Institute of Technology

Abstract : This paper describes the chaotic properties observed in the 2-mass shear system incorporating
elasto-plastic force-displacement relation. If a small amount is added to the acceleration at a certain time inthe system, the value of variation of the subsequent acceleration may diﬀer signiﬁcantly. The orbital expansion
rate of the system was calculated, and it was found that the orbital expansion rate was distributed with fourdelta peaks. We ﬁnd that the orbital expansion rate of one of the delta peaks is positive, which allows us toobserve chaotic properties. Furthermore, when the term of damping in this system is set to 0, the distributionof the orbital expansion rate becomes a distribution with one peak in the positive value, indicating that theterm of damping has the eﬀect of suppressing the appearance of chaotic properties.
Keywords : chaotic properties, 2-mass shear system, elasto-plastic force-displacement relation, time history
response analysis

1. はじめにため，実際の建物の振動との関係だけを考えるならば建物の地震応答解析によく利用されるモデルに 2 質

図 2 の方がイメージしやすいかもしれない．このモデ点系せん断モデルがある．各層の質量を質点に集中さルについての運動方程式を立て，外力として地震動をせ，層せん断力がその層の変形だけにしか効かないと入れることで質点の運動が決まり，それにより建物のすることで，層の剛性を 1 つの独立なバネで表し，各各層の振動状態を知ることができる．

質点を直列につないだ図 1 のようなモデルである．図このモデルは，地震動が小さい場合は各層の振動状

1 を見ると，縦振動を連想してしまうかもしれないが，せん断モデルなので質点は水平方向に振動する．その態をよく表すが，地震動が大きくなると，実際の建物の振動状態とのずれが大きくなる．地震動が大きいと

第一工業大学研究報告 第32号 （2020）

k2
m2

m2
k2
k1

m1

m1

動加速度を ẍ0〔m/s2 〕としたとき()() ()()
m1 0ẍ1
c11 c12ẋ1+0 m2ẍ2
c13 c14ẋ2)()()() (ẍ0k11 k12x1
m1 0+=0 m2ẍ0k13 k14x2k1(1)が解くべき運動方程式である．ただし，剛性マトリク

図 1: 2 質点系せん断モ

図 2: 2 質点系せん断モ

デル

デルのイメージ

各層の振動が大きくなり，建物が変形し始めるためである．建物の変形が大きくなると，きれつ，降伏，すべりなどの現象が生じ，復元力と変形の関係は履歴ループを描くようになる（弾塑性復元力特性）．このようスは，各バネの剛性 ki〔N/m〕を用いて)) ((k1 + k2 −k2
k11 k12=−k2k2
k13 k14となり，また，減衰マトリクスは初期剛性比例型で与えるため(な状態になると，このモデルから得られる各層の振動

c11c12状態は現実と大きくずれ始める．しかしながら，弾塑

c13c14

性復元力特性をモデルに組み込むことで，これが改善される．弾塑性復元力特性をモデルに組み込む方法は多々提案されており，弾塑性復元力特性を 2 質点系せん断モデルに組み込むことで，変形が大きくなっても各層の状態を知ることができる．
本論文は，地震応答解析を行う上で重要となる 2 質)2h=ω(k1 + k2−k2−k2k2)(3)となる．ただし，h は減衰定数であり，ω は 1 次固有円振動数である．
弾塑性復元力特性については，図 3 のようなバイリニア型復元力特性と図 4 のようなスリップ型復元力特性を 1:1 の割合で組み込む．

Q点系せん断モデルに弾塑性復元力特性を組み込んだ運動方程式自体に着目する．2 質点系せん断モデルの運γk動方程式に，弾塑性復元力特性を組み込むと，それまで観測できなかった興味深い性質が観測されるようになる．その性質とは，ある時刻にわずかな摂動を加えることで，その後の数値計算結果が全く違う値となるkOxδという性質である．これはカオスの性質である可能性が高い．
本論文では，2 質点系せん断モデルに弾塑性復元力2δ特性として，バイリニア型復元力特性とスリップ型復元力特性を 1 : 1 で組み込み，そこに外力として地震動を入れた運動方程式のもつカオス的1 性質について報告

図 3: バイリニア型復元力特性

する．

Q
2. 扱う運動方程式
本論文では，2 質点系せん断モデルを扱う．2 質点系せん断モデルは，図 1 のように，質点 2 つが縦につながり，それぞれの質点が水平方向に振動するモデルだが，実際に解くべき運動方程式は以下のようになる．質点の質量を mi〔kg〕（ただし，i = 1，2 であり，i = 1 が 1O階を表し，i = 2 が 2 階を表す．
），位置を xi〔m〕，地1 本研究は，外力として地震動を運動方程式に入れる．地震動は持続時間が有限のため，その時間内でカオスと断定するのは難しい．
そこで，有限時間で得られる数値計算結果でカオスかどうかを判定
し，「カオス的」と表現する．

(2)

図 4: スリップ型復元力特性x山田・福永：２質点系せん断モデルに弾塑性復元力特性を組み込んだ系で観測されるカオス的性質伏変位 δi〔m〕は同じにしている．また，図 3，図 4 の縦軸は弾塑性復元力 Q〔N〕，横軸は位置 x〔m〕である．この弾塑性復元力特性により，式 (1) の数値計算を行う際，各時刻の質点の状態に合わせて剛性マトリクスの値が変化する．これにより，建物に変形が生じて加速度 ẍ0〔cm/s2 〕

なお，バイリニア型復元力特性とスリップ型復元力
特性の初期剛性 ki〔N/m〕，第 2 剛性 γki〔N/m〕，降80
40200−20−40−60−80も各層の状態をよく表すことができる．

基礎部分

60

本論文では，数値積分法の 1 つである Newmark の020406080

100120

時間 t〔s〕

β 法（β = 0.25）を用いて式 (1) を数値計算する．
図 5: 基礎部分で計測された加速度

3. 数値計算結果と実データとの比較
まずは，式 (1) の数値計算が実際の建物の振動状態8 時 48 分ごろ発生した地震（震源地：日向灘，マグニチュード：6.3）時に，鹿児島県霧島市に実在する木造 2階建てで計測された加速度データと数値計算結果の比較を行う．この建物には，基礎部分および 2 階天井部分に加速度計が設置してあり，地震時の加速度が記録されている．この加速度データを用いて数値計算を行う．基礎部分の加速度データを地震波の加速度として式 (1) に入力し，2 階天井部分で計測された加速度と質点 2 の加速度（数値計算）の比較を行う．数値計算時

加速度〔cm/s2 〕をよく表すことを確認するために，2019 年 5 月 10 日

80

2 階天井
数値計算

6040200−20−40−60−80020406080100120

時間 t〔s〕
図 6: 2 階天井で計測された加速度と数値計算結果に必要となる入力データは表 1 の通りである．表 1 中の剛性比については，初期剛性に対する第 2 剛性比である．また，重力加速度 g = 9.80665 m/s2 ，減衰定数4. 摂動を加えたときの加速度の変動
式 (1) が実際の建物の振動状態をよく表すことが確認できたところで，式 (1) の数値計算に対して，ある

表 1: 入力データ

時刻の加速度に摂動を加える．数値計算アルゴリズム

1階

2階やパラメータ，入力データはすべてそのままの状態で，各階重量 Mi〔kgf〕

10000

7000

ある時刻の加速度 a〔cm/s2 〕に微小量 ∆a〔cm/s2 〕を初期剛性 Ki〔kN/rad〕

2000

1500加えて，その後の加速度の変動を数値計算する．

降伏変形角 γyi 〔rad〕

0.010.01

図 7 は，図 6 の数値計算結果（青線）および数値計算途

剛性比 γki0.015

0.015

中の時刻 t = 40.0 s のときに，加速度 a〔cm/s2 〕に微小量 ∆a = 0.092 cm/s2 を加えて計算した結果（赤線）で

h = 0.13，数値計算時の刻み幅 ∆t = 0.01 s とした．刻ある．なお，微小量 ∆a〔cm/s2 〕は，−0.1 5 ∆a 5 0.1

み幅 ∆t については，計測された加速度のサンプリングの一様乱数として与えた．∆a〔cm/s2 〕を加えた直後は間隔と同じにした．なお，数値計算に必要な降伏変位差が見えないが，t = 70.0 s を過ぎた辺りから差が見え始3
δ = γyi × H〔m〕と初期剛性の単位変換 ki = Ki ×10Hめる．図 8 は，図 7 の t = 39.8 s から t = 42.0 s を切り出〔N/m〕で必要な仮想階高 H は 3.0 m とした．

した図である．t = 40.0 s のときに，∆a = 0.092 cm/s2

図 5 が実際に計測された基礎部分の加速度である．縦が加えられていることが見て取れる．その後，1 秒弱程軸が地震波の加速度と仮定している基礎部分の加速度度は微小な差が見えるが，その後は差が無くなってい2ẍ0〔cm/s 〕，横軸が時間 t〔s〕である．また，図 6 が2 階天井部分で計測された加速度と数値計算の結果をした図である．図 9 を見ると，t = 40.0 s に ∆a〔cm/s2 〕重ねたグラフである．図 6 を見ると，数値計算の結果を加えた数値計算結果と加えていない数値計算結果がと実際に計測されたデータは完全に一致するわけでは明らかに違う．t = 40.0 s で ∆a〔cm/s2 〕を加え，そのる．図 9 は，図 7 の t = 72.2 s から t = 85.0 s を切り出ないが，P 波部分（t = 0 から 40 秒過ぎまでの加速度後，1 秒弱程度で値の差が無くなり，その状態が 30 秒が小さい部分）の加速度，S 波部分（P 波以降の部分）程度続いた後，再び差が表れ始める．しかも再び現れの加速度ともによく一致している．

た差は，図 8 の t = 40.0 s から t = 40.5 s でみられる微第一工業大学研究報告 第32号 （2020）

80

加速度〔cm/s2 〕

加速度〔cm/s2 〕

80

数値計算
∆a を加えて数値計算

6040200−20−40−60−80

数値計算∆a を加えて
数値計算

6040200-20-40-60

020406080

100-8012002040時間 t〔s〕

数値計算
∆a を加えて数値計算

0.40.2
0−0.2
−0.4

40.5

4141.50.40.2
0-0.2

4240

40.5

100−10−20788042

図 11: 弾塑性復元力特性なしの系の t = 39.8〜42.0

加速度〔cm/s2 〕

2076

41.540

数値計算∆a を加えて数値計算

7441時間 t〔s〕

図 8: ∆a を加算したときの t = 39.8〜42.0

加速度〔cm/s2 〕

120

数値計算
∆a を加えて数値計算

0.6

時間 t〔s〕

30100-0.4404080

図 10: 弾塑性復元力特性なしの系に ∆a を加算

加速度〔cm/s2 〕

加速度〔cm/s2 〕

図 7: t = 40.0 で ∆a を加算

0.660

時間 t〔s〕

8284

数値計算∆a を加えて数値計算

3020
100-10-207274

時間 t〔s〕

7678808284

時間 t〔s〕

図 9: ∆a を加算したときの t = 72.2〜85.0

図 12: 弾塑性復元力特性なしの系の t = 72.2〜85.0

小な差ではなく，明らかに違う値を取っている2 ．

を見ると，図 9 と違い完全に一致している．つまり，あこのような結果を引き起こす原因を特定するために，る時刻に微小量を加算したときに，図 9 のような，全弾塑性復元力特性を組み込まずに数値計算を行ってみく違う変動を引き起こす原因は弾塑性復元力特性の影る．図 10 は，弾塑性復元力特性を組み込まずに数値計響によるものだとわかる．

算を行った結果である．先ほどと同様に t = 40.0 s のときに微小量 ∆a〔cm/s2 〕を加えて計算した．図 10 を見ると，∆a〔cm/s2 〕を加えても加えなくても結果が一致していることがわかる．図 11 は，図 10 の t = 39.8 sから t = 42.0 s を切り出した図である．図 8 と同様，

t = 40.0 s で ∆a〔cm/s2 〕が加算されているのが見て取れる．その後，1 秒程度は差の部分が所々で見られるが，その後は一致している．また，図 12 は，図 10 のt = 72.2 s から t = 85.0 s を切り出した図である．図 122 この結果は，加える微小量により出現しないこともある．

5. 軌道拡大率
弾塑性復元力特性を運動方程式に組み込むことで，数値計算途中の加速度 a〔cm/s2 〕に微小量 ∆a〔cm/s2 〕を加えると，その後，全く違う変動をすることがある．
微小変化が急速に拡大し，全く違う挙動を示す変動はカオスの性質である可能性が高い．カオスかどうかを判断するためには，まずリアプノフ指数を計算し，その正負を調べるのが一般的である．しかしながら，地震動を外力として入力している運動方程式においては，

山田・福永：２質点系せん断モデルに弾塑性復元力特性を組み込んだ系で観測されるカオス的性質

地震波が有限のため，軌道拡大率の極限を取ることが0.1できず，リアプノフ指数を計算することができない．そ動をする部分において軌道拡大率を計算する．
軌道拡大率は，微小量 ∆a が n ステップ後にどれだ

0.08

相対度数

こで，微小量 ∆a〔cm/s2 〕を加えたときに全く違う変け大きくなるかの指標であり，最初の差を ∆a0 とし次0.06
0.040.02の差を ∆a1 としたときに∆a1= f ′ (a0 )∆a00-0.05

-0.04

-0.03

(4)-0.02

-0.01

00.01

軌道拡大率 Λと f ′ を定め，さらに次のステップで ∆a1 が ∆a2 になっ

図 13: 軌道拡大率の分布たとすると∆a2= f ′ (a1 )∆a1∆a2= f ′ (a1 )f ′ (a0 )∆a0(5)(6)考えられる減衰の項を 0 にしたときの軌道拡大率 Λ を計算する．減衰の項を 0 にするためには，式 (3) の減衰定数 h = 0 とすればよい．先ほどと同様，t = 40.0 sのとき，−0.1 5 ∆a 5 0.1 の一様乱数を加え，アンサとなり，これを n ステップ繰り返すとn−1∏∆anf ′ (xt )=∆a0t=0

最後に，カオス的性質を妨げる要因になっているとンブル数 105 で軌道拡大率 Λ を計算する．計算の結果，(7)

軌道拡大率 Λ = 0.0041 と正の値が得られた．また，図14 は，h = 0 での軌道拡大率 Λ の分布である．図 13 ととなる．この右辺を0.1′nΛf (xt ) = e(8)t=0とおくことで軌道拡大率 Λ を定義する．式 (8) を見ると，Λ > 0 のときは，n → ∞ で微小量が指数関数的に0.08

相対度数n−1∏

0.06
0.040.02発散することがわかる．また，Λ < 0 のときは，微小量が指数関数的に小さくなり，n → ∞ で 0 となる．な0-0.05

お，具体的に軌道拡大率 Λ を数値計算するときは，次-0.04-0.03-0.0100.01

軌道拡大率 Λの式を使う．

1∑ln |f ′ (xt )| >n t=0

図 14: 減衰定数 h = 0 での軌道拡大率の分布n−1<Λ> = <-0.02(9)は違い，Λ = 0.0035 にピークを 1 つもつような分布が

ただし，< · · · > はアンサンブル平均を意味する．この得られた．図 15 は，減衰定数 h = 0 で，t = 40.0 s の式 (9) の n → ∞ の極限を取ったときの収束値 Λ∞ がリときに ∆a〔cm/s2 〕を加えたときの加速度 a〔cm/s2 〕アプノフ指数であり，Λ∞ > 0 はカオスの性質の 1 つの変動である．図 7 のときより，早い段階で値の差がに求め，その値 Λ でカオス的性質をもつかを判断する．
軌道拡大率 Λ の計算は次のように行った．t = 40.0 sのときに，−0.1 5 ∆a 5 0.1 の一様乱数を加え，アンサンブル数 105 で軌道拡大率 Λ を計算した．軌道拡大率Λ を計算した結果，軌道拡大率 Λ = −0.0044 と，Λ < 0という結果が得られた．そこで，105 個の軌道拡大率 Λの分布を取った．図 13 は軌道拡大率 Λ の分布である．
図 13 を見ると，Λ = 0.0050，− 0.0025，− 0.013，− 0.046のところにデルタピークをもつような分布が得られた．

4 つのデルタピークのうち 1 つが Λ = 0.0050 > 0 と正の値を取る．このとき，系はカオス的性質をもち，図9 のような全く違う変動が観測される．

加速度〔cm/s2 〕である．ここでは，有限時間の軌道拡大率 Λ を数値的
100806040200-20-40-60-80020406080100120

時間 t〔s〕
図 15: 減衰定数 h = 0 での加速度 a の変動大きくなった．図 16 は，図 15 の t = 80.0〜100.0 s の第一工業大学研究報告 第32号 （2020）

ときを切り出した図であるが，全く違う変動をしていカオス的性質が観測されるが，それ以外ではカオス的る様子が見てわかる．

性質は観測されない．与える微小量と Λ の関係を明確にすることは，カオス的性質の出現条件の解明につな

加速度〔cm/s2 〕

3010がる．また，減衰定数 h = 0 から h を増加させるとき，Λ の分布がどのように変化していくのか，2 つ目，3 つ目，4 つ目のデルタピークがどのように発生するのか，0そのメカニズムを解明していく必要もある．

20-10●-20-3080859095100

時間 t〔s〕

カオス的性質が数値解析に及ぼす影響について

2 質点系せん断モデルに弾塑性復元力特性を組み込んだ運動方程式は，カオス的性質をもつために，数値解析を行う際は注意が必要である．系がカオス的性質をもつということは，十分小さな差が，その後の変動に

図 16: 減衰定数 h = 0 での加速度 a の変動（t = 80.0

大きな影響を与えるということである．実際，カオス的〜100.0）性質が出現しないような状況では，微小量 ∆a〔cm/s2 〕の影響は無視できたが，そうでない状況では，加速度

6. 考察

a〔cm/s2 〕に微小量 ∆a〔cm/s2 〕を加えると，その後の変動が大きく異なっていた．これは，数値積分を行

カオス的性質の出現についてう際，十分気をつけなければならないことである．例

2 質点系せん断モデルにおいて，剛性に弾塑性復元力えば，どの数値積分アルゴリズムを採用するかにより，特性を組み込まなかった場合，カオス的性質が現れな結果が大きく異なってくる．また，どの地震波を利用●かったことから（図 10，11，12 参照），カオス的性質するかでも解析結果が変わる．また，同じ地震で計測が現れるのは剛性に弾塑性復元力特性を組み込んだとされたデータを利用する場合でも，計測器の性能の違きであることがわかった．また，図 13，14 より減衰のいによる誤差で解析結果が大きく変わってしまう．数項がカオス的性質の出現を抑える効果があることがわ

値解析を行う際は，系がもつカオス的性質を十分考慮かった．このことから，図 7 において，t = 40.0 s のとしながら解析を行う必要がある．

きに微小量 ∆a を加えても，相対的に加速度 a〔cm/s2 〕
塑性復元力特性を組み込んだ系でも変位が小さいため

7. まとめ
本論文は，2 質点系せん断モデルに弾塑性復元力特初期剛性のみの運動となる）カオス的性質は現れない．

性を組み込んだ系で観測されるカオス的性質について

その後，加速度 a〔cm/s2 〕の値が相対的に大きくなる記述した．2 質点系せん断モデルは，弾塑性復元力特性S 波部分になると，弾塑性復元力特性の効果により，カオス的性質が出現する．図 8 において，t = 40.0 s でを組み込むことによりカオス的性質が観測されること加えられた微小量 ∆a〔cm/s2 〕は，1 秒程度経つと無効果があることも明らかになった．の値の小さい P 波部分では，建物の変形がないため（弾がわかった．また，減衰の項がカオス的性質を抑えるくなったように見えるが，実は非常に小さい値として今後は，軌道拡大率 Λ と与える微小量の関係，減衰残っており，その小さな値が S 波部分でカオス的性質定数 h を変化させたときの Λ の分布の変化について研により大きな値となって現れることで，図 9 のような究を進めていく．

全く違う変動が観測される．
このカオス的性質は実際の建物にも含まれる性質であると考えられ，条件によってはカオス的性質が現れると考えられる．実際の建物においてカオス的性質が現れる条件は何か，どのような影響を及ぼすかを明確

謝辞
本研究を遂行するにあたり，貴重なデータを提供していただいた西日本工業大学 デザイン学部 建築学科
古田智基 博士，横浜国立大学 大学院都市イノベーショ
ン研究院 特別研究教員 中尾方人 博士に深く感謝する．にすることで，耐震，制振，免振の技術向上に役立つと考えられる．
●

軌道拡大率の分布について
図 13 を見ると，軌道拡大率 Λ のデルタピークが 4 つ

参考文献[1] 柴田明徳,“最新耐震構造解析（第 3 版）”, 森北出版
株式会社, 2014

ある．これは，与える微小量 ∆a〔cm/s2 〕により Λ が[2] 井上政義, 秦浩起，“カオス科学の基礎と展開 −複どの値を取るかが決まることを表す．Λ > 0 であれば雑系の理解に向けて−”, 共立出版株式会社, 1999