計算工学講演会論文集 Vol. 28 (2023 年 5 月)

計算工学会

粒子法と ASI-Gauss 法のカップリングに基づく
波浪–骨組構造連成解析
Wave–Framed Structure Interaction Analysis Based on
Coupled Particle Method and ASI-Gauss Method
大村浩之,1) 三目直登,2) 浅井光輝,3) 磯部大吾郎 4)
Hiroyuki Ohmura, Naoto Mitsume, Mitsuteru Asai and Daigoro Isobe
1)2)

博 (工) 筑波大学 システム情報工学研究群 (〒 305-8573 茨城県つくば市天王台 1-1-1, E-mail: homura@bosai.go.jp)

博 (工) 筑波大学 システム情報系 助教 (〒 305-8573 茨城県つくば市天王台 1-1-1, E-mail: mitsume@kz.tsukuba.ac.jp)
3)4)

博 (工) 九州大学 工学研究員 准教授 (〒 819-0395 福岡県福岡市西区元岡 744, E-mail: asai@doc.kyushu-u.ac.jp)

博 (工) 筑波大学 システム情報系 教授 (〒 305-8573 茨城県つくば市天王台 1-1-1, E-mail: isobe@kz.tsukuba.ac.jp)

In this reseach, we developed a partitioned FSI analysis method by coupling both stabilized ISPH method
and ASI-Gauss method aimed for tidal wave and framed structures interaction problems. For transmissionof physical quantities between ﬂuid and structure domains, we adopted a particle-beam element interaction
model based on Improved ERP wall boundary model. Some test calculations were conducted as veriﬁcationand validation for the proposed method. Furthermore, we showed results of simulation that simple framed
structures were broken and swept away by a ﬂow mimicking tsunami.
Key Words : ISPH method, FEM, ASI-Gauss method, Fluid–structure interaction, Tsunami

1. 緒言
津波による構造物の被害を予測するためには，津波によって構造物が変形，破壊する現象を把握することが肝要である．そのためには数値解析の利用が不可欠だが，既存の研究例としては，津波遡上解析 [1] や，構造物に衝突する波浪の解析 [2] など，流体のみの挙動を再現する解析への適用に留まっているのが現状である．
すなわち，津波によって構造物が流出するようなダイ
ナミックな連成問題を一貫して再現する手法は確立されていない．そのような手法が実現すれば，津波が家屋等を巻き込みながら遡上する過程など，より複雑で現実に近い状況を再現でき，津波対策のさらなる高度化が期待できる．
上記の背景を踏まえ，本研究では粒子法とはり要素の有限要素法（FEM）をカップリングさせた分離型解法の Fluid-Structure Interaction (FSI) 解析手法を構築し，波浪によって骨組構造物が変形，破壊，流出する現象
を一気通貫でシミュレーションする．
粒子法には非圧縮性流れに対して定式化された Incompressible Smoothed Particle Hydrodynamics (ISPH) 法 [3]を採用する．粒子法では点群で連続体を離散化するため，境界における幾何学的な情報が失われる．そのため，一般に境界条件を付与するために何らかのモデリングが必要となる．壁面上の境界条件処理に関しては，
連成面上での物理量計算に直結するため，壁境界モデルの精度は連成解析に適用する上で特に重要である．そこで本研究では，高精度に壁面境界条件を処理でき，かつ連成計算の精度が優れる Improved Explicit Represented
Polygon (IERP) モデルを開発し，これを ISPH 法に適用した．

構造解析に適用するはり要素の FEM としては，Adaptively Shifted Integration (ASI)-Gauss 法 [4] に基づく解
析スキームを用いる．ASI-Gauss 法は線形 Timoshenko
はり要素における数値積分点と応力評価点の位置関係[5] を利用し，弾塑性変形の計算を高精度化した手法である．さらに，本研究では構造物が流出し壁面および他の構造物に衝突する様子を再現するために，梁部材の部材形状を厳密に考慮した接触モデル，ならびに破断面を陽に表現する破断モデルを提案する．
粒子法とはり要素のカップリングとしては，まず梁部材の形状を補間するポリゴンを定義し，それを粒子法における壁境界として扱う．その上で，節点とポリゴン頂点，ならびに要素とポリゴン面それぞれで評価される物理量を関連づけることで，粒子–はり要素間の連成モデルを構築する．

2. ISPH 法に基づく自由表面流れ解析(1) 離散化スキーム
本研究では次式で記述される非圧縮性流れを対象と
する．
∇·v=0
(1)
Dv
1
= − ∇p + ν∇2 v + g(2)
Dtρ
ただし，v は流速，t は時間，ρ は密度， p は圧力，ν は動粘性係数， g は重力加速度である．流体解析スキー
ムには ISPH 法 [3] を用いる．ISPH 法をはじめとする
SPH 系列の手法では，次式のカーネル近似に基づく粒子離散表現を基本とする．
X mj
ϕ(xi ) ≈ ⟨ϕ⟩i =ϕ j Wi j
(3)
ρj
j∈P
i

ここで，ϕ は任意の物理量，i, j は粒子のインデックス，
x は位置ベクトル，mi , ρi は粒子 i 上で定義される質量と密度，Wi j は粒子 i, j 間の距離に依存するカーネル関数である．また，Pi は粒子 i の影響半径 h 内に存在する粒子の集合である．カーネル関数には 3 次スプライン関数を用いる．
導出過程は割愛するが，式 (3) を用いて次式の 1 階，
2 階の空間微分モデルが導かれる．
X mj ∇ϕ(xi ) ≈ ⟨∇ϕ⟩i =ϕ j − ϕi ∇Wi j
(4)
ρj
j∈P
i

X mj 
 ri j∇2 ϕ(xi ) ≈ ∇2 ϕ = 2
ϕ j − ϕi· ∇Wi jiρj|ri j |2
j∈P
D

E

(5)

i

ここで， ri j = x j − xi である．
非圧縮性流れを対象とする粒子法では，一般に射影法 [6] に基づいて時間方向の離散化を行う．すなわち，
速度と圧力を変数分離し，速度は陽的に，圧力は陰的に計算する．n + 1 ステップ目の圧力 p は次式のポアソン方程式を解くことで得る [7]．

∇2 pn+1 =α (ρ − ⟨ρ⟩n )ρ
∇ · v∗ +
∆t
∆t2

(6)

ここで，v∗ は中間状態における流速，∆t は時間増分，ρ̂は離散化された状態での数値的な流体密度である．また，α は安定化係数で，時間増分の値に応じて設定される．

(2) IERP 壁境界モデル
IERP モデルおよびその前身である ERP モデル [8] では，図–1 に示すように，壁面をポリゴンで近似的に表現し，粒子 i をポリゴン k に関して鏡映したゴースト粒
子 ik が流体粒子から受ける寄与を計算し，それをもう一度鏡映することでポリゴン k が粒子 i に与える壁寄与を計算する．例えば 1 階微分量の壁寄与 ⟨∇ϕ⟩Si,k は次式のように求められる．
X mj 
ϕ jk − ϕi ∇Wik j (7)⟨∇ϕ⟩Si,k = Ri,k ⟨∇ϕ⟩Fik = Ri,kρjj∈P

角が凸形状（π ≤ θk ≤ 2π）の場合，





∀l ∈ Wnear
Wik jk ; rl j · nl < 0Ŵik j = 

0
else

(9)

となる．ここで，θk は角の角度，Wnear
は角 k に隣接すkるポリゴンの集合， rl j はポリゴン l から粒子 j に向かう最短経路方向ベクトル， nl はポリゴン l 上の外向き
法線ベクトルである．また，凸形状の角に隣接するポリゴンに関しても同様の考えに基づいてカーネル関数が修正される．
以上に記した方法によって，物理量の壁寄与項が表現される．例えば圧力勾配項は次式のように与えられる．
* +S X 
!






pi
∇p

 X m j p jk
k
=
+
R
∇
Ŵ
(10)

i,k
i j




ρ i
ρ
ρ
ρ


j
j
i
k∈Wi j∈P ki

ここで，Wi は粒子 i の影響半径内に存在するポリゴンと角の集合である．同様に，粘性項は
!
* +S X 






pi
∇p
 X m j p jk

k
=
+
R
∇
Ŵ
(11)

i,k
i j




ρ i
ρ
ρ
ρ


j
j
i
k∈Wi j∈P kiとなる．3 次元空間への拡張など，IERP モデルのより詳細な説明は文献 [9] を参照されたい．

図–1 ポリゴン面に関する鏡映に基づく間接的な壁寄与表現

ik

ここで，ik は粒子 i をポリゴン k に関して鏡映した粒子である．Ri,k はその鏡映を表す表現行列である．詳細は割愛するが，0, 2 階微分項の寄与についても同様に計算される．このように壁寄与を計算することにより，ポリゴンで表現された特定の壁領域が粒子に与える力を陽に記述できる．すなわち，その反作用を壁面に与えれば境界面上の運動量保存が保証されるため，既存の壁境界モデルと比べて連成解析の精度および安定性が向上する．
さらに，ポリゴン同士が接する角に対しても同様に寄与を考慮することで，より精緻に壁領域を表現する．
角 k が粒子 i に与える寄与もポリゴンの場合とまったく同様に計算されるが，カーネル関数 Wik j は壁領域の体積を整合させるために次式に従って修正される．角が凹形状（0 ≤ θk < π）の場合，

Ŵik j =2π − 3θk
Wi k jθk(8)

3. ASI-Gauss 法に基づく骨組構造解析(1) 離散化スキーム
空間方向の離散化として，全ての変位成分を同一の 1
次 Lagrange 関数で補間する線形 Timoshenko はり要素
[10] を用いる．また，時間積分には Newmark のβ法を用いる．これらを用いて時空間的に仮想仕事の原理式
を離散化すれば，次式の変位 dt+∆t を未知数とした連立方程式が得られる．
M d̈t+∆t + C ḋt+∆t + Kdt+∆t = f t+∆t − gt
(12)!
1  t+∆t−dt 
1 t
1 t
d̈t+∆t =d−d̂ + 1 −
d̈
(13)
2
β∆t
2β
β∆
!
!
γ t
γ
γ  t+∆t−dt 
t+∆t
+ 1−
d
d̂ + 1 −∆t d̈t (14)
ḋ
=
β
2β
β∆2
ここで， M, C, K はそれぞれ系全体の質量マトリック
ス，減衰マトリックス，剛性マトリックスである．また， f , g はそれぞれ節点外力ベクトル，節点内力ベク
トルで，γ, β は Newmark のβ法の重みパラメータである．本研究では，γ = 5/6, β = 4/9 とする．
ASI-Gauss 法では，梁部材を 2 つの線形 Timoshenko
はり要素で表現したうえで，弾性領域ではそれぞれの応力評価点を 3 次はり要素のものと一致するように配置
し，塑性ヒンジが生じた後はその位置に応力評価点を配置するように積分点をシフトする．これにより，通常のFEM に比べて弾塑性変位解の精度が著しく向上する．

(2) 梁部材接触モデル
梁部材同士および梁部材と壁面との接触現象を表現するために，本研究では粒子ベース接触モデル（図–2）を提案する．以下では，当該モデルにおける接触力計
算の流れを記す．
まず，実際の部材の断面形状を考慮して，はり要素 iの部材表面を ΓSi と定義する．いま，はり要素 i に接触する要素を j とすると，幾何学的には線分として表されるはり要素 i, j 間での最近傍点は 1 点に定まる．その最近傍点を中心とし，部材の幾何形状に応じて規定される接触が起きうる領域を接触領域 Γci とする．
部材形状を厳密に考慮して接触力を計算するために，
接触領域 Γci 上に等間隔に接触粒子群 Ci を配置し，粒子同士の接触力計算を介して節点上での接触力を与えnる．接触粒子 a ∈ Ci が b ∈ C j から受ける衝突力 fab
は，
力積法 [11] に基づいて次式のように計算される．
 (1+e)v ·n msab ab ab(vab · nab < 0)nab + fab−n∆t
fab = 
(15)

fs
else
ab

ただし，e は反発係数，vab , mab は粒子 a, b 間で計算さsは数値誤差等に起因れる相対速度および相対質量， fab
する粒子貫入を防ぐための安定化力である．また，接触面の単位法線ベクトル nab は，na − nbnab ≡|na − nb |(16)

na は粒子 a ∈ Ci 上での ΓSi の外向き単位法線ベクトルである．
t同様に，摩擦力 fab
は，!
(1 + e)vab · t ab mabtnfab= min µ fab
, −
t ab
(17)
∆t
ただし，µ は摩擦係数である．また，摩擦方向に対する接触面の単位接線ベクトル t ab は，t ab ≡ −vab − (vab · nab ) nab|vab − (vab · nab ) nab |

(18)

である．
c
n
t
接触力 fab= fab+ fab
によって生じる要素軸回りのトc
ルク mab は次式のように求められる．
noccmcab = (xa − xi,a ) × fab− ( fab· ei )ei(19)
ここで， xa , xi,a はそれぞれ粒子 a の絶対位置および要
素にマッピングした位置，ei は要素軸方向の単位ベクトルである．

以上のように計算した接触力および接触トルクは次式に従って要素 i の節点力ベクトル fic に変換される．
!
c
1 XX ⊤fabcfi =N (ξi,a ) c
(20)
mab
Nab
a∈Ci b∈C j

ただし，Nab は総接触点数，N は形状関数マトリックス，
ξi,a は xi,a を無次元化した量である．総接触点数で除しているのは，多点接触の簡易的なモデリングで，全ての接触点の寄与が同一であることを仮定している．

Element

Contact
domainContact
particle
Element

図–2 粒子ベース接触モデル(3) 梁部材破断モデル
本研究では，構造物が津波によって流出する挙動を追跡するために，梁部材の破断面を陽に表現する陽的破断モデルを開発し導入する．
はじめに，梁部材の破断を扱った既存研究 [12] を参考に，次式のひずみ量を参照した判定式を用いて破断の発生判定を行う．
!2
!
!2
!2
!
χy 2
εyz 2
χx
εzz
ε xz
+
+
+
+
≥1
(21)
χcx
χcy
εczz
εcxz
εcyz
ここで，χ x , χy , εzz , ε xz , εyz はそれぞれ要素座標系 x, y 軸方向の曲率，軸ひずみ，x, y 軸方向のせん断ひずみである．また，上添え字の c はひずみの破断臨界値を表す．
はり要素 i が破断した時刻を t f ，その時の応力評価
点および数値積分点の無次元位置をそれぞれ r f , s f とする．破断面では，破断が発生した瞬間に要素が保持している断面力 σt f (r f ) が解放力 g f として作用する．両者の関係は，
g f = B⊤ (s f )σt f (r f )(22)となる．ただし， B は節点量を要素量に変換するマトリックスである．
破断面を陽に表現するために，破断した瞬間に r f で要素を分離させ，節点を生成する．同時に，破断面に接する節点に式 (22) の解放力を作用させる．その際，破断過程におけるエネルギー減衰を模擬して解放力を平滑化する．平滑化には，次式で表される破断時刻 t f からの経過時間 ∆t f に依存する近似 Heaviside 関数 H(∆t f )
を用いる．
(
!)
∆t f
1
H(∆t f ) =1 + cos π(23)2τf
ここで，τ f は平滑化パラメータで，本研究では経験的に τ f = ∆t f /100 とする．

4. ISPH 法と ASI-Gauss 法のカップリングに基づく波浪–骨組構造連成解析
波浪–骨組構造連成解析を実現するために，上述した
ISPH 法と ASI-Gauss 法をカップリングした分離型解法の FSI 解析手法を構築する．時間発展アルゴリズムには Conventional Serial Staggered (CSS) アルゴリズム[13] を採用する．以下では粒子とはり要素の間における運動量伝達の計算方法，すなわち連成モデルについて記述する．
(1) 構造から流体に対する連成モデル
はり要素の節点上で評価される物理量から，連成境
界面，すなわち部材表面上の量を計算する．IERP モデ
ルでは部材表面を含めた壁面をポリゴンで表現するため，ポリゴンの頂点での量を求めればよい．
構造から流体に対する連成モデルの模式図を図–3 に示す．部材表面を表すポリゴンの頂点 M の位置 x M と，
頂点 M に対応する節点 B の位置 xB の間には次の関係
式が成り立つ．


t+∆t
xt+∆t
= Pt+∆txtM − xtB
(24)
M − xBB
PB は節点 B における断面の回転に対する表現行列で，節点の回転方向の変位に基づいて決定される．
頂点 M の速度および加速度については，変位 ut+∆t
M =
t+∆t
t
x M − x M に Newmark のβ法を適用することで計算できる．ただし，その際の重みパラメータには構造解析と同じ値を用いる．

素軸回りのトルク mFi,k→S は，


 n

 o
beam
F →S
F →S
mFi,k→S = xwall− fi,k· ee ee (26)i,k − xi,k,e × fi,k
ポリゴン k に対応するはり要素 e に作用する節点力
feF →S は次式のように書ける．
F →S !
XXfi,k
F →S⊤ wallfe=
N (ξi,k ) F →S
(27)
mi,k
i
k∈Wi

wall
ただし，ξi,kは xbeami,k,e を無次元化した要素座標である．
上述した流体から構造に対する連成モデルにおける計算は，ポリゴン k が 2 次元角もしくは 3 次元角に置き換わっても同様に成り立つ．

Front view

Top view

Particle:

Polygon:Element:

図–4 流体から構造に対する連成モデル(3) 妥当性検証
以上で提案した連成解析手法の精度・妥当性検証として，梁部材の水中落下問題やゴム板付きダムブレイク問題を解く．紙面の都合上，計算結果の詳細は発表にて述べるが，良好な精度で流体構造連成現象を再現できることが確認された．

図–3 構造から流体に対する連成モデル(2) 流体から構造に対する連成モデル
流体から構造に対する連成モデルの模式図を図–4 に示す．IERP モデルを用いて壁領域を表現すると，粒子i がポリゴン k に及ぼす力は次式のように計算できる．
 * +SDE  ∇pF →S
2 S
= −mi −fi,k+ ν ∇ v 
(25)
i,k
ρ i,k
F →S
流体力 fi,kの作用点は粒子 i から見たポリゴン k 上
wallの最近傍点 xi,k とする．また，xwall
i,k をはり要素 e 上に
F →S
投影した座標を xbeamとする．
fによって生じる要i,k,ei,k5. 波浪および漂流物の衝突を受ける建物の流出解析
本章では，波浪–骨組構造連成問題に対する提案手法の有用性を示すために，津波を模した段波およびそれに伴う漂流物の衝突を受ける建物が変形・破壊し，流出する過程をシミュレートした結果を示す．
(1) 計算条件
図–5 に解析対象とする系の概観を示す．実際の津波と建物のスケールを想定し，十分な幅方向長さ，奥行4 [m]，高さ 8 [m] の閉空間の中に，幅 5 [m]，奥行 4 [m]，高さ 1.2 [m] の水を配置する．また，水の端から 3 [m]
離れた位置に 2 層 1 × 1 スパンの単純骨組構造物（以下，建物 A と呼称する）を配置する．さらに，建物 A
から 3 [m] に建物 A と同一の構造物（以下，建物 B と呼称する）を配置する．建物 A, B のスパン長および階
高はそれぞれ 2 [m] とし，建物を構成する全ての部材の断面形状は 0.24 [m] 角の正方形とする．

Top viewY

Unit: [m]
1
2

4
1
X
Front view

2
8
2

5
Z

3

1.2

3
Building A

2
Building B

X
図–5 系の概観
表–1 流体解析条件

Time increment
Number of time steps
Particle spacing
Number of particles
Eﬀective radius
Relaxation coeﬃcient
Density
Kinetic viscosity

5.0 × 10−4 [s]10000
5.0 × 10−2 [m]
192000
1.2 × 10−1 [m]5.0 × 10−3
1.0 × 103 [kg/m3 ]1.0 × 10−6 [m2 /s]

表–2 構造解析条件

Time increment
Number of time steps
Young’s modulus
Poisson’s ratio
Yield stressDensity
Diameter of contact particle
Coeﬃcient of restitution
Coeﬃcient of friction
Smoothing parameter5.0 × 10−4 [s]10000
7.93 × 107 [N/m2 ]0.4
5.98 × 105 [N/m2 ]4.91 × 102 [kg/m3 ]2.0 × 10−2 [m]
0.15
0.3
−2
1.0 × 10 [s]

段波は予め入力した速度時刻歴に従って左端の壁面を動かすことで発生させる．0.0 ≤ t < 0.1 [s] では線形的に速度を増加させ，0.1 ≤ t < 0.9 [s] の区間では一定，
0.9 ≤ t < 1.0 [s] では線形的に減少させ，1.0 [s] 以降は常に速度を零とする．また，最大速度は 3 [m/s] とする．
流体解析条件および構造解析条件をそれぞれ表–1 と
表–2 に示す．構造物の全ての部材は 2 要素で分割している．構造物の材料はスギ材を想定して物性値を決定している．また，ここでは木材の脆性的な破壊特性を考慮し，部材が塑性化した瞬間に破断も生じるものとして，式 (21) の破断臨界値は降伏ひずみと同値とする．

(2) 計算結果
図–6 に建物 B のみを配置した場合（以下，ケース 1と呼称する）および建物 A, B 両方を配置した場合（以下，ケース 2 と呼称する）それぞれの可視化結果を示
す．流体粒子のコンターは圧力 [Pa]，はり要素のコン
ターは降伏関数を表している．
ケース 1 では，建物 B の 1 階部分の柱が段波の衝突を受けて変形を起こすものの，塑性化・破壊には至らなかった．一方，ケース 2 では，最初に段波の衝突を受ける建物 A の 1 階における上流側の柱に大きな流体力が作用し，底面および 2 階との接合部分で破断が生じた．その後，波は後流側の 1 階柱にも衝突し，底面との接合部で破断を引き起こした．1 階の柱を全て失った
建物 A は転倒し，その衝撃によってほとんどの接合部で破断が生じた．漂流物と化した建物 A は流れに乗って移流し，崩壊直後に建物 B に衝突した．衝突を受けた建物 B の 1 階における上流側の柱は底面との接合部で破断を起こし，建物 B も建物 A と同様のプロセスで崩壊を起こした．
以上の計算結果から，津波および建物のスケールにおいて，波浪の衝突を受けた骨組構造物が変形・破壊を起こし，流出することで他の構造物に衝突する過程を提案手法によって再現できることが示された．また，建物単体では被害が生じない場合においても，上流側で流出した構造物が漂流物として衝突を起こすことによって，被害を伝播させる可能性があることが示唆さ
れた．

6. 結言
本研究では，粒子法における壁境界モデルとして連成解析に適した IERP モデル，ならびに梁部材の流出と衝突を再現できる粒子ベース接触モデルと陽的破断モ
デルをそれぞれ開発した．その上で，それらを導入した ISPH 法と ASI-Gauss 法をカップリングした波浪–骨組構造連成解析手法を構築した．計算例として，波浪および漂流物の衝突を受ける建物の流出解析を実施した．計算結果から，漂流物の衝突を介した被害伝播過程を含む波浪による骨組構造物の被害シミュレータとして提案手法が有用であることが示された．
参考文献[1] Kowalik, Z. and Murty, T.S.: Numerical simulationof two-dimensional tsunami runup, Marine Geodesy,
Vol.16, No.2, pp.87–100, 1993.
[2] Kleefsman, K.M.T. et al.: A volume-of-ﬂuid based
simulation method for wave impact problems, J. ofComp. Phys., Vol.206, No.1, pp.363–393, 2005.
[3] Cummins, S.J. and Rudman, M.: An SPH projection
method, J. of Comp. Phys., Vol.152, No.2, pp.584–
607, 1999.
[4] Lynn, K.M. and Isobe, D.: Finite element code for impact collapse problems of framed structures, Int. J. forNum. Methods in Eng., Vol.69, No.12, pp.2538–2563,
2007.
[5] Toi, Y.: Shifted integration technique in onedimensional plastic collapse analysis using linear andCase 1Case 2

図–6 可視化結果（左：ケース 1，右：ケース 2）

cubic ﬁnite elements, Int. J. for Num. Methods in Eng.,
Vol.31, No.8, pp.1537–1552, 1991.
[6] Koshizuka, S. and Oka, Y.: Moving-particle semiimplicit method for fragmentation of incompressible
ﬂuid, Nuc. Sci. and Eng., Vol.123, No.3, pp.421–434,
1996.
[7] Asai, M. et al.: A stabilized incompressible SPH
method by relaxing the density invariance condition,
J. of Applied Maths., Vol.2012, pp.1–24, 2012.
[8] Mitsume, N.: Explicitly represented polygon wallboundary model for the explicit MPS method, Comp.
Particle Mechs., Vol.2, No.1, pp.73–89, 2015.
[9] 大村浩之: 粒子法と ASI-Gauss 法のカップリングに基づく流体–骨組構造連成解析, 博士論文, 筑波大学,
2023.
[10] Hughes, J.R.T. et al.: A simple and eﬃcient ﬁnite element for plate bending, Num. Methods in Eng., Vol.11,
No.10, pp.1529–1543, 1977.

[11] Baraﬀ, D.: Interactive simulation of solid rigid bodies, IEEE Comp. Graphs. and Apps., Vol.15, No.3,
pp.63–75, 1995.
[12] Isobe, D. and Tanaka, S.: Sequential simulationsof steel frame buildings under multi-phase hazardous
loads during earthquake and tsunami, Frontiers in
Built Env., Vol.7, pp.669601, 2021.
[13] Farhat, C. and Lesoinne, M.: Two eﬃcient staggered
algorithms for the serial and parallel solution of threedimensional nonlinear transient aeroelastic problems,Comp. Methods in Applied Mechs. and Eng., Vol.182,
No.3-4, pp.499–515, 2000.