計算工学講演会論文集

Vol.29 (2024年6月)

日本計算工学会

地震動を受ける木造家屋の倒壊シミュレーションCollapse Simulation of a Wooden House under Seismic Excitation伊藤慎也1)，大村浩之2) ，磯部大吾郎3)
Shinya Ito, Hiroyuki Omura and Daigoro Isobe
1) 筑波大学大学院（〒305-8573 茨城県つくば市天王台1-1-1，Email: s2302834@u.tsukuba.ac.jp）
2) 博（工）防災科学技術研究所 特別研究員（〒305-0006 茨城県つくば市天王台3-2, E-mail: homura@bousai.go.jp）
3) 博（工）筑波大学 システム情報系 教授（〒305-8587 茨城県つくば市天王台1-1-1, Email:isobe@kz.tsukuba.ac.jp）
We simulated the collapse of a wooden house using the ASI-Gauss method. A method to reduce thestiffness of wood beams based on the plastic region of the element cross section was used as a constitutive
law of wood beams. The validity of the numerical method was verified by comparing with a shaking table
experiment. In addition, the seismic excitation observed during the 2024 Noto Peninsula Earthquake wasused as input for the analysis.
Key Words : ASI-Gauss method, Constitutive law of wood beams, Wooden house, Collapse simulation

1． 序論と一致するように数値積分点を配置するという特徴を有木造家屋の倒壊を含む非線形地震応答を解くためのシしており，これを用いると最小分割数の2要素分割でも収ミュレーション手法として，T. Nakagawa and M. Ohta,束解が得られる．さらに，部材の破断および部材間での接2003[1][2][3]の Extended Distinct Element Method（EDEM）

触判定が容易であるという特長を有する．

をベースとした手法が広く知られている．EDEMでは部材をバネ・ダッシュポットで接続された剛体球の集合で表

3． 木材構成則現するため，部材の破壊に伴う不連続な変位場の扱いに(1) 木材はりの特徴は適しているが，連続体理論の枠組みで部材の弾塑性挙考慮すべき木材の特徴として，圧縮力による塑性ひず動を扱うことができない．

みに対し，引張力による塑性ひずみが小さいことが挙げ

そこで本研究では，骨組構造の崩壊挙動を含む強非線られる．長谷川ら[7]は，木材はりの弾塑性解析において形構造解析に有効なASI-Gauss法[4]にはり部材の木材構Bilinear型の応力-ひずみ曲線を用いることには妥当性が成則を導入する．岩崎ら[5]によって示された1方向の曲げあるとしている．そこで本研究では，図-1に示すBilinear型と軸力を考慮した木材構成則を拡張し，2方向の曲げと軸の応力-ひずみ曲線を用いる．図中の𝜎 は降伏応力，𝜀 は力を考慮可能な木材構成則を提案する．これにより木造引張破壊ひずみ，𝜀 は圧縮破壊ひずみである．

家屋などの三次元空間での弾塑性解析に木材構成則が導(2) 正規化入できるようになり，塑性域の広がりによる剛性低下と本研究では，はりの断面を長方形断面とする．さらに，破壊を表現できるようになる．提案した木材構成則を応力分布を求めるために，要素断面を1×1の正方形とな

ASI-Gauss法に導入することで，連続体としての木材の弾るように正規化する．これに伴い断面力を以下の式で正塑性変形から破壊までの過程と，破壊後の飛散といった規化する．

離散物としての過程を統一的にシミュレートすることが可能となる．この提案手法を用い，2005年に実大三次元振動破壊実験施設（以下，E-ディフェンス）で実施された既𝑚⬚𝑀, 𝑚⬚
𝜎 𝑏ℎ𝑀,𝑛𝜎𝑏 ℎ𝑁𝜎 𝑏ℎ(1)存木造家屋倒壊実験[6]の再現解析を行うことで，提案手法の妥当性と有用性を示す．また，能登半島地震で記録された地震動を入力波形として解析を行う．

2． ASI-Gauss法[4]本研究で用いたASI-Gauss法は，線形チモシェンコはり要素を用いて弾塑性解析を行う有限要素法（FEM）の一種

これにより，正規化平面では降伏応力の値は-1となる．また，曲げモーメントの正負によって応力分布が鏡像になるので，議論を簡単にするするため曲げモーメントは常に正であると仮定する．正規化された断面の弾性域の応力分布を式(2)で定義する．

である．弾塑性解析時には数値積分点を順応的にシフトすることで，少ない要素分割数で高精度な解析が可能となる．特に，弾性域での応力評価点が部材のガウス積分点

𝑓 𝑥, 𝑦𝑎𝑥𝑏𝑦𝑐 𝑎0, 𝑏0(2)𝑭𝟎となる𝑎, 𝑏, 𝑐が得られる．このように，重積分で得られる断面力と構造解析から得られた断面力が一致する
𝑎, 𝑏, 𝑐を求め，式(3)に代入することで弾塑性境界𝑔 𝑥 を求める．

図-1 応力-ひずみ曲線

図-2 応力分布と弾塑性境界

𝐹 𝑎, 𝑏, 𝑐𝑚

𝑚 𝑎, 𝑏, 𝑐(6-1)

𝐹 𝑎, 𝑏, 𝑐

𝑚𝑚 𝑎, 𝑏, 𝑐

(6-2)

𝐹 𝑎, 𝑏, 𝑐𝑛⬚𝑛⬚ 𝑎, 𝑏, 𝑐

(6-3)

𝑎𝑏𝑐

𝑱𝑎𝑏𝑐

𝑱

𝜕𝐹⎡⎢ 𝜕𝑎
⎢𝜕𝐹
⎢ 𝜕𝑎
⎢𝜕𝐹
⎣ 𝜕𝑎

𝜕𝐹
𝜕𝑏
𝜕𝐹
𝜕𝑏
𝜕𝐹
𝜕𝑏

𝐹 𝑎, 𝑏, 𝑐
𝐹 𝑎, 𝑏, 𝑐
𝐹 𝑎, 𝑏, 𝑐𝜕𝐹⎤𝜕𝑐 ⎥𝜕𝐹 ⎥𝜕𝑐 ⎥𝜕𝐹 ⎥𝜕𝑐 ⎦(7-1)(7-2)

図-3 塑性パターン(4) 断面積・断面二次モーメントの算出
Bilinear型の応力-ひずみ曲線を用いたため，塑性域では1とな る． 弾塑性 境界 𝑔 𝑥 は ，式(2)をメントを算出する．断面二次モーメントについては，初め1とすることで式(3)のように求められる．応力に図心を求め，図心を基準とした断面二次モーメントを

常に 𝑓 𝑥, 𝑦𝑓 𝑥, 𝑦

求められた弾塑性境界をもとに断面積・断面二次モー算出する．この計算式も，塑性パターンごとに異なる．

分布と弾塑性境界を図-2に示す．

(5) 引張・圧縮破壊判定
𝑎𝑥𝑏𝑔 𝑥𝑐

弾性域の応力分布をもとにひずみ分布を求める．図-2の

1(3)

𝑏

断面で最も圧縮ひずみが大きい位置は(0,0)であり，引張ひずみが大きくなる位置は(1,1)である．(0,0)のひずみが圧縮破壊ひずみ𝜀 になった際に圧縮破壊，(1,1)のひずみが(3) 応力分布の算出
応力分布は断面力から直接求めることができないが，断面力は式(4)に示すように応力分布を重積分することで容易に求められる．その計算結果は，図-3に示す塑性パターンによって異なる．例として塑性パターン１の場合の重積分の結果を式(5)に示す．

𝑚 𝑎, 𝑏, 𝑐

文献[5]で行われた解析と同じのモデル（図-4，表-1）で𝑓 𝑥, 𝑦 ∗ 𝑦(4-1)

𝑓 𝑥, 𝑦 ∗ 𝑥

1𝑑𝑆
2(4-2)

を𝑥軸周りにそれぞれ0,15,30,45°回転させて解析を行う．

(4-3)[mm]，支持はりの高さは100 [mm]である．

2方向の曲げモーメントを作用させるため，中央はり部分

𝑓 𝑥, 𝑦 𝑑𝑆⬚

𝑚 a, b, c𝑛 𝑎, 𝑏, 𝑐

2𝑎 𝑏
3𝑎 𝑏

解析を行う．コの字のはりの上部に水平荷重Pを加える．

なお，はりの断面形状は100×100 [mm]，スパンは1000

𝑛 𝑎, 𝑏, 𝑐
𝑚 𝑎, 𝑏, 𝑐

4． 要素解析

1𝑑𝑆
2⬚2𝑎𝑏行う．

(1) 解析モデル⬚𝑚 𝑎, 𝑏, 𝑐

引張破壊ひずみ𝜀 になった際に引張破壊を判定する．ただし，後述する再現実験では別の方法を用いて破断判定を

2𝑏 1 𝑐24𝑎𝑏1

2𝑎 1 𝑐24𝑎 𝑏3𝑎𝑏16𝑎𝑏1𝑐𝑐𝑐6𝑎𝑏𝑐(2) 解析条件

(5-1)

要素実験では木材構成則に基づいて塑性化を表現する
ため，数値積分点のシフトは行わない．

(5-2)(3) 解析結果
はり中央の𝑧軸方向のたわみを図-5に示す．解析結果か(5-3)

ら，角度をつけるほど塑性化開始時の変位が小さくなることがわかる．これは，図心から圧縮端までの距離が長く，Newton-Raphson法に用いる関数𝑭を式(6)のように設定し，式(7)のように逐次的に𝑎, 𝑏, 𝑐を求めていくことで，圧縮端でのひずみが大きくなるためである．
塑性化開始後の剛性に注目する．角度をつけているモデルは塑性化開始が早いが塑性化後の剛性低下が緩やか

図-4 解析モデル
表-1 木材の材料パラメータ
𝜎 : 𝑐𝑜𝑚𝑝𝑟𝑒𝑠𝑠𝑖𝑣𝑒 𝑠𝑡𝑟𝑒𝑛𝑔ℎ 𝑀𝑃𝑎

47.5𝜀 : 𝑡𝑒𝑛𝑠𝑖𝑙𝑒 𝑓𝑎𝑖𝑙𝑢𝑟𝑒 𝑠𝑡𝑟𝑎𝑖𝑛 108.56

𝜀 : 𝑐𝑜𝑚𝑝𝑟𝑒𝑠𝑠𝑖𝑣𝑒 𝑓𝑎𝑖𝑙𝑡𝑢𝑟𝑒 𝑠𝑡𝑟𝑎𝑖𝑛 101.49

𝐸: 𝑦𝑜𝑢𝑛𝑔 𝑠𝑚𝑜𝑑𝑢𝑙𝑢𝑠 𝐺𝑃𝑎

9.560

図-5 荷重-𝑧軸方向変位いて破断判定を行う．
𝜃𝜃∨ｃ𝜃

∨𝜎𝜎(8)である．これは，角度のあるモデルでは塑性域が断面の角から進行し，同じ曲率でも欠損部分が比較的小さくなる

ここで，𝜃 , 𝜃 はそれぞれ要素座標系𝑥, 𝑦軸方向の変形角，からである．また，角度をつけるほど破壊判定に至る荷重

𝜃 , 𝜃 は変形角の破断臨界値，𝜎, 𝜎 はそれぞれ軸方向応が小さいことが確認できる．これは，角度をつけるほど圧力と引張強度である．終局変形に達した際の靭性破壊お縮・引張端までの距離が大きくなり，それに伴って端部でよび引張強度に達した際の脆性破壊の2通りの破壊パタのひずみが大きくなるからである．これらの結果から，塑

ーンを考慮する．

性域の広がりによる剛性低下が妥当に再現できているこ(2) 解析モデルとが確認された．

木造家屋の解析モデルは，全ての部材を2要素分割のはり要素で作成する．本研究では簡単のため壁材と柱を一5． 移築した既存木造住宅の倒壊実験

体の構造要素としてモデル化し，構造要素試験で得られ2005年にE-ディフェンスで実施された木造家屋倒壊実た荷重変位関係と一致するように断面二次モーメントと験[6]では，実際に建っていた既存不適格の木造家屋を震塑性断面係数をキャリブレーションした．1階および2階動台上に移築し，JR鷹取波100％による加振が実施された．の梁および屋根に該当する要素の密度は，解析モデルの実験では同等の仕様を持つ2棟の家屋を移築し，片方を無重量が5節で述べた重量の実測値と一致するように決定補強のまま，もう片方を耐震補強したうえで加振実験がしている．弾性自由振動解析で算出された解析モデルの1実施されたが，本稿では無補強試験体のみを解析対象と次固有周期はX（EW）方向が0.347 [s]，Y（NS）方向が0.523する．

[s]である．X，Y方向ともに実測値に対して値が大きいの

実験の過程で実測された対象家屋の重量は屋根60.62

は，外壁および内壁の初期剛性を過小評価しているため[kN]，2階133.90 [kN]，1階101.29 [kN]，総重量295.82 [kN]と考えられるが，地震動下では早期に剛性低下が生じるであった．また，1次固有周期はX（EW）方向が0.137 [s]，ため，倒壊挙動への影響は少ないと考えられる．

Y（NS）方向が0.222 [s]であった．

要素の材料定数を表2に示す．断面寸法は柱が10×10[cm]，梁が10×15 [cm]もしくは10×20 [cm]，筋交いが 3

6． 木造家屋倒壊実験の再現解析(1) 復元力特性モデル
本研究では，木材の復元力特性を図-6に示すようにモデ×9 [cm]，その他の部材は全て10×15 [cm]である．
(3) 解析条件
時間増分は0.001 [s]，Newmark 𝛽法におけるパラメータ4/9，𝛾5/6とする．物理減衰は導入していない．

ル化する．図中の弾塑性領域では，3節で述べた木材構成は𝛽則モデルを用い，断面内の塑性（欠損）部分を同定し剛性接触は部材と床面の間のみ考慮する．解析の入力波形にを低下させる．要素解析では，応力が最大強度に達したらはJR鷹取波100％加振時に震動台上で実測された加速度破壊状態と判定していたが，実際の架構は一定の耐力を時刻歴を用いる．

保持するため，再現解析では，完全塑性領域で塑性流れ則(4) 解析結果に従うと仮定した．除荷時は常に初期剛性で挙動するも図-7に再現解析の可視化結果および実験映像を示す．コのとする．この復元力特性モデルは，一般的なTri-linearモンターは圧縮強度で正規化した最大圧縮応力を表していデルの第2直線領域を木材構成則に基づく非線形領域にる．また，破断し耐力が喪失した要素は可視化していない．

置き換えたものといえる．

図から，4.4 [s]で家屋がY軸負方向へ大きく傾き，6.0 [s]に終局状態は要素の破断で表現するものとし，式(8)を用は1階の部分的な損傷が生じ，7.6 [s]の段階には1階の層崩図-6 木材の復元力特性モデル
表-2 木材の材料パラメータ
Hinoki

Larch

Westernhemlock

Young’s modulus [GPa]

8.8

9.8

10.3

Poisson’s ratio

0.4

0.4

0.4

Compressive strength [MPa]

39.0

44.0

40.0

Tensile strength [MPa]

117.0

132.0

120.0

Maximum strength [MPa]

85.0

85.0

85.0

0.16

0.16

0.16

440500480

Ultimate deformation angle[rad]Density [𝑘𝑔/𝑚 ]図-7 解析の可視化結果（左）と実験映像（右）壊が完了している．実験結果と比較すると，層崩壊が始めるまでの過程は概ね一致しているが，実験では層崩壊の進行速度が解析よりも遅い．これは，耐力を失った部材が図-8 1 階における Y 軸方向の層間変位時刻歴

2階部分に対して緩衝材のような効果を発揮しているためと考えられる．
1階におけるY軸方向の層間変位時刻歴を図-8に示す．
グラフより，倒壊挙動の定性的な傾向が実験結果と対応していることが確認できる．

7． 能登半島地震を入力波形とした解析(1) 解析モデル・解析条件
解析モデルは6節と同じモデルを用いる．解析の入力波形には能登半島地震で記録されたK-NET富来波100％の加速度時刻歴を用いる．その他解析条件は6節と同じである．
(2) 解析結果
図-9に可視化結果を示す．図から17.2 [s]では14.9 [s]と図-9 解析の可視化結果

最大圧縮荷重がほぼ同じにも関わらず，1階が大きく変形している様子が確認できる．これは，16 [s]付近で塑性化を起こし，その後除荷されたが塑性ひずみが残ったことが原因だと考えられる．図-10に示した1階におけるY軸方向の相関変位時刻歴からも同じことが推測できる．17.8[s]で1階に部分的な損傷が生じ，21.2 [s]には家屋が完全に崩壊している．

8． 結論
本研究では既存の研究を拡張し，2方向の曲げモーメン

図-10 1 階における Y 軸方向の層間変位時刻歴

トと軸力が作用する木材はりの弾塑性解析手法を開発しJ. Wood Sci. 49, 499–504, 2003.

た．要素解析の結果から，塑性域の広がりによる剛性低下[3] T. Nakagawa, M. Ohta, T. Tsuchimoto and N. Kawai,と破断を妥当に表現できることが確認された．そして木Collapsing process simulations of timber structures under

材構成則を導入したASI-Gauss法に基づく木造家屋倒壊

dynamic loading III: numerical simulations of real-size

シミュレーション技術を開発し，E-ディフェンスで実施

wooden houses. J. Wood Sci. 56, 284–292, 2010.

された既存木造家屋倒壊実験の再現解析を実施し妥当性[4] 磯部大吾郎：はり要素で解く構造動力学：建物の崩壊を検証した．また能登半島地震を入力波とした解析を行

解析からロボット機構の制御まで，丸善出版，2020.
[5] 岩崎正二，佐々木寛幸，山本亮，宮本裕：塑性域の広った．

がりを考慮した平面木骨組構造物の弾塑性解析，第参考文献

３回木橋技術に関するシンポジウム論文報告書，95–[1] T. Nakagawa and M. Ohta, Collapsing process simulationsof timber structures under dynamic loading I: simulationsof two-story frame models. J. Wood Sci. 49, 392–397, 2003.
[2] T. Nakagawa and M. Ohta, Collapsing process simulations
oftimberstructures

underdynamicloadingII:

simplification and quantification of the calculating method.

102, 2004.
[6] 槌本敬大, 移築した既存木造住宅の耐震補強・無補強実験. 建築技術 675, 158–161, 2006.
[7] 長谷川薫，薄木征三：曲げとせん断を受けるひずみ軟化はりの弾塑性変位解析，土木学会論文集，351–358, 1989.