計算工学講演会論文集 Vol. 27 (2022 年 6 月)

計算工学会

ISPH 法と ASI-Gauss 法のカップリングに基づく
流体構造連成解析手法の開発
Development of Fluid-Structure Interaction Analysis Method
Based on Coupling ISPH Method and ASI-Gauss Method
大村浩之 1) 三目直登 2) 浅井光輝 3) 磯部大吾郎 4)
Hiroyuki Omura, Naoto Mitsume, Mitsuteru Asai and Daigoro Isobe
1)修 (工) 筑波大学 システム情報工学研究群（〒 305-8573 茨城県つくば市天王台 1-1-1, E-mail: s2030186@s.tsukuba.ac.jp）
2)

博 (工) 筑波大学 システム情報系 助教（〒 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 research, we developed a partitioned FSI analysis method with a coupling stabilized ISPH method
and ASI-Gauss method aimed for tsunami and buildings interaction problems. For transmission of physicalquantities between ﬂuid and structure domains, we adopted a particle-beam element interaction model basedon Improved Explicitly Represented Polygon (IERP) wall boundary model. According to the numericalresults simulating a free falling square column member into a water obtained by the developed method, it
is conﬁrmed that the balance between buoyancy and gravity, and rotation angle of the ﬂoating member at
equilibrium state can be reproduced.
Key Words : Incompressible SPH Method, Polygon Wall Boundary Model, ASI-Gauss Method, Linear
Timoshenko Beam Element, Fluid-Structure Interaction, Partitioned Coupling Scheme

1. 緒言
津波による沿岸地域の被害を軽減するために，津波避難ビルや津波ハザードマップの整備が進められている．
津波避難ビルの設計指針 [1] では，津波外力を予想される波高に準じた静水圧の 1.5〜3 倍に定めている．すな
わち，動的な波圧成分を単純化して設計外力を定めている．また，漂流物の衝突など，より複雑な荷重は考慮されていないのが現状である．津波ハザードマップ
については，浅水長波方程式による津波遡上解析に基づいていることが多い．その際，建物などの構造物は粗度係数などの簡便な形で表現もしくは無視されて地形がモデリングされるため，津波によって建物が破壊され，瓦礫が流れに巻き込まれたりするなどの複雑な現象は再現できない．
上述の複雑な現象を考慮し，津波避難ビルやハザー
ドマップの信頼性を向上させるためには，流体構造連成（Fluid-Structure Interaction, FSI）シミュレーションの利用が不可欠である．既存の FSI 解析手法は数多く存在するが，津波のような自由表面流れを含む FSI 問題に好適であるとして，粒子法ベースの FSI 解析手法
[2][3] が注目されている．これらは流体と構造を全て粒子法で解く一体型解法の手法 [2] と，流体は粒子法で解
き，構造解析は有限要素法（FEM）など別のスキームで解く分離型解法の手法 [3] に類別される．ところが，構造物を粒子や一般的な有限要素で離散化する場合，建物のような大規模構造物に対しては計算コストの観点から適用性が低い．したがって，実スケールで津波と建物の連成問題を解くための FSI 解析手法が望まれる．

また，近年は「津波減災」の概念が重要視されており，減災設計の観点では漂流物の衝突を含む津波外力下での終局限界荷重や崩壊挙動などの把握が必要と思われる．つまり，津波-建物連成問題に適用される FSI の構
造解析スキームには，崩壊挙動を含む強非線形解析に対応できることが求められる．
以上を踏まえ，本研究では，津波-建物連成問題を対象とした，粒子法とはり要素を用いた FEM のカップリングに基づく分離型解法の FSI 解析手法を開発する．粒子法には安定化 ISPH 法 [4]，構造解析には ASI-Gauss
法 [5] を用いる．安定化 ISPH 法は半陰解法の粒子法である ISPH 法の圧力ポアソン方程式に安定化項を導入したスキームで，圧力解の精度・安定性と体積保存性を両立できる．ASI-Gauss 法は線形 Timoshenko はり要素を用いた FEM に基づく骨組構造解析手法で，数値積分
点をシフトすることにより，はり要素内の任意の位置で塑性ヒンジや破断面の表現が可能な手法である．すなわち，少ない要素分割数で高精度な塑性崩壊荷重解が得られるため，大規模骨組構造物の崩壊解析に適している．
また，開発する分離型解法 FSI 解析手法には，Improved Explicitly Represented Polygon (IERP) 壁境界モデル [6] をベースとした粒子-はり要素連成モデルを適用する．IERP 壁境界モデルは，分離型解法の FSI における連成モデルに適した壁境界モデルとして知られている ERP モデル [7] を半陰解法に適用できるよう高精度化したモデルである．
本稿では，安定化 ISPH 法および ASI-Gauss 法の基礎理論と粒子-はり要素連成モデルについて説明した後，開発する連成解析手法の計算フローについて述べる．最後に，開発手法を用いていくつかのベンチマーク問題
を解き，その妥当性を検証する．

2. 解析手法(1) 安定化 ISPH 法に基づく非圧縮自由表面流れ解析a) 支配方程式と時間積分スキーム
ラグランジュ記述の非圧縮流れの支配方程式は次式のように書ける．
Dv
1= − ∇p + ν∇2 v + g
Dt
ρ

(1)

ρ は密度，v は流速，p は圧力， g は重力加速度，ν は動粘性係数である．半陰解法の粒子法では，射影法に基づき，時間ステップ間の中間状態を定義して速度を段階的に更新する [8]．中間状態における仮速度 v∗ は，次式のように粘性項と重力加速度の作用によって決まる．
v∗ − vn= ν∇2 vn + g∆t(2)

n は時間ステップ，∆t は時間増分である．仮速度をさらに圧力勾配項の作用によって更新することで，新しいステップでの流速を得る．
vn+1 − v∗
1= − ∇pn+1
∆t
ρ

(3)

安定化 ISPH 法 [4] では，次式に示す圧力ポアソン方
程式を解くことで圧力を得る．

∇2 pn+1 = −ρ0ρ0 − ρn∇ · v∗ + α
∆t
∆t2

(4)

右辺第二項は粒子分布の乱れを防ぐように付加的な圧力を与える安定化項で，α は緩和係数である．α の値は∆t1/3 に設定することが適切であると数学的に示されている [9]．
b) SPH 法による離散化
SPH 法 [10] は，積分形式のカーネル近似式を出発点とし，それを有限個の計算点（粒子）を用いて重み付
き総和の形で次式のように表す．
∑ mj
ϕ(xi ) ≈ ⟨ϕi ⟩ =W(xi j , h)ϕ j
(5)
ρj
j∈P
i

⟨⟩ は離散化された量であることを表す．ϕ は任意の物
理量，i, j は粒子番号，m, ρ, x はそれぞれ粒子の質量，密度，位置ベクトル，W, h はそれぞれカーネル関数とカーネルのサポート半径（影響半径）である．Pi は粒子 i の影響半径内に存在する粒子の集合である．なお，
ϕi j = ϕ j − ϕi としている．本稿では，カーネル関数として 5 次スプライン関数を用いる．
式 (5) を用いると，式 (1)(4) における各項の離散化式が得られる．本研究では以下の離散化モデルを使用
する．
 





⟨∇pi ⟩ ∑   p j pi m j  2 + 2  ∇Wi j 
(6)
=





ρi
ρ
ρ
j∈Pi

j

i

⟨ν∇2 vi ⟩ =∑ ( ρ j ν j + ρi νi xi j · ∇Wi j )mjvi jρ j ρi|xi j |2j∈P⟨∇2 pi ⟩ =

∑(

mj

j∈P

⟨∇ · vi ⟩ =ρ j + ρi xi j · ∇Wi jpi j
ρ j ρi|xi j |2

∑ (mj
j∈P

ρj

(7)

)

(8)

)
vi j · ∇Wi j(9)c) その他の解析技術
自由表面粒子の検出には Tanaka and Matsunaga[11] の近傍粒子数に基づく方法を用い，検出しきい値は 0.80とする．また，DS 法 [12] を導入し，粒子間距離が一定の距離（初期粒子間距離の 90 ％）を下回った粒子ペアには，その距離に応じた力積を与えて粒子分布を補正する．
(2) ASI-Gauss 法に基づく骨組構造解析a) 支配方程式と時間積分スキーム
有限要素法で離散化された連続体の運動方程式を次式に示す．
Mü + Cu̇ + Ku = F − FR
(10)
M, C, K はそれぞれ質量マトリックス，減衰マトリック
ス，剛性マトリックスである．ü, u̇, u はそれぞれ節点加速度ベクトル，節点速度ベクトル，節点変位ベクトルである．F, FR はそれぞれ節点外力ベクトルと節点残差
力ベクトルである．式 (10) を Newmark’s β 法で時間方向に離散化すると，次式を得る．
(
)
1
δ
M
+
C
+
K
un+1 =
β∆t
β∆t2
{(
)
}
1
1 n
Fn+1 − FnR + M− 1 ün +
u̇
(11)
2β
β∆t
{(
)
(
) }
δ
δ
+C
− 1 ∆tün +− 1 u̇n2ββ
ここで，β, δ は解の安定性や精度に影響するパラメータである．本研究では，計算の安定性を向上させるために，数値減衰作用が得られる β = 4/9，δ = 5/6 の設定[13] を採用している．
b) 線形 Timoshenko はり要素を用いた離散化
本研究では建物のような骨組構造物の解析を効率よく行うために，はり要素で構造物を離散化する．はり要素には Timoshenko 梁理論に従う線形 Timoshenko はり要素（Linear Timoshenko Beam Element, LTBE）[14]を用いる．LTBE が表現する断面力は二軸の曲げモーメ
ント，軸力，ねじりモーメント，二軸のせん断力で，変位は曲率，ねじり角，三軸の並進方向変位の計 6 自由度である．変位ベクトルの要素内分布は線形変位関数
によって次式のように表される．

u(ξ) =

1
1(1 − ξ)u1 + (1 + ξ)u2
2
2

(12)

u1 , u2 はそれぞれ要素のローカル節点番号 1, 2 における節点変位ベクトルである．ξ は要素上の正規化座標
で，−1 ≤ ξ ≤ 1 の値をとる．詳細な導出は省略するが，

Timohenko 梁理論に式 (12) を適用することで，剛性マ
トリックス K などが得られる．
質量マトリックス M ははり要素の分布質量マトリッ
クスとして求められる．定式化の詳細は文献 [5] などで説明されている．
減衰マトリックス C は，レイリー減衰に基づいて次式のように計算する．
C = αM + βK

Fluid domain

Polygon(13)

比例係数 α, β は物理減衰特性を決定するモデルパラメー
タである．

Wall domainc) 数値積分点のシフティング
Toi[15] は LBTE と剛体ばねモデル [16] のエネルギー
等価性から，次式を導いた．

Fig.1 Schematic of a wall domain expressed by mirror particles

s1 = −r1

ただし，Ri,k は粒子 i に対して定義される壁 k に関する
鏡映写像の表現行列（以下，鏡映行列）で，次式のように計算できる．

(14)

上式は，数値積分点の位置を変えることで任意の位置の応力を評価できることを意味する．ASI 法 [17] では，この関係を利用して，降伏が検出された位置に応じて
数値積分点をシフトすることで，要素分割を細かくせずとも骨組構造物の塑性崩壊荷重を精度よく得られる．
また，ASI-Gauss 法 [5] では 1 部材を 2 つの LTBE でサブセットとして表現し，サブセット全体に対して 2 点
積分を適用するようにそれぞれの要素の数値積分点を配置することで，弾性変位解を 3 次の精度で解くことができる．数値積分点のシフトを含む計算フローの詳
細は文献 [5] に記載されている．

3. 粒子-はり要素連成モデル(1) IERP 壁境界モデル
粒子とポリゴン壁の相互作用モデリングには，Improved Explicitly Represented Polygon (IERP) 壁境界モデル [6] を用いる．IERP 壁境界モデルでは，注目する流体粒子 i が壁から受ける寄与 ϕwall
i,k を，粒子 i を壁に′関して鏡映した粒子 i を用いて表現する．壁として扱う対象は，ポリゴンで表される平面と，2 つのポリゴンが接する線分として定義される 2 次元角，3 つ以上のポリゴンが交わる点として定義される 3 次元角である．
なお，IERP モデル自体の V ＆ V は文献 [6] で既に行っている．
a) 平面に対する鏡映′
Fig.1 に示すように，ミラー粒子 j でポリゴン壁 k の内側領域が表現されているとする．このとき，例えば粒子 i に作用する圧力勾配項の壁面寄与成分は式 (6) から次のように書ける．







⟨∇pi ⟩wall,k ∑ pi 
 ′  p j′
(15)=
m j  2 + 2  ∇Wi j′ 





ρi
ρ
ρ
′
′
j ∈Pi

j

i

ここで，任意の粒子の組み合わせ i, j に対して，以下の等式が成り立つxi j′ = Ri,k xi′ j(16)∇Wi j′ = Ri,k ∇Wi′ j(17)

Ri,k = I − ni,k ⊗ ni,k

(18)

ni,k は平面ポリゴン k の粒子 i に向かう単位法線ベクト′′ルである．また，m j = m j , ρ j = ρ j などが成り立つ．したがって，式 (15) は次式のように書き換えられる．
 





⟨∇pi ⟩wall,k ∑ pi   p j′
m j  2 + 2  Ri,k ∇Wi′ j 
=





ρi
ρ
ρ
j
i
j∈P ′

(19)

i

鏡映粒子上の圧力 p j′ は，壁面 k 上の圧力ノイマン条件
を満たすように決定される．

{(
}
)
wall
p j′ = p j + 2ρ j |x j − xnear· n j,k n j,kj,k | g − ai,k(20)

ここで，xnearj,k はポリゴン k 上の粒子 j に対する最近傍点
wall
位置，a j,k は xnear
j,k におけるポリゴン k の加速度である．
式 (19) から分かるように，IERP モデルでは粒子 i に′対する壁面寄与成分はその鏡映粒子 i と近傍の流体粒
子 j ∈ Pi′ の情報のみで記述できる．すなわち，仮想的な粒子を生成することなく，壁面境界条件を厳密に満たすような壁面寄与成分を与えることができる．
式 (7–9) の他の離散化式についても，同様の方法で壁′面寄与成分を定式化できる．その際，鏡映粒子 j の速度は，壁面 k の速度ディリクレ条件に基づいて次式のように計算できる．
(
)

(
)

R j,k v j + 2 vwall
free slip
j,k · n j,k n j,k(21)
v j′ = ()
 − v j + 2vwall
no slipj,k角に対する鏡映IERP モデルでは，壁領域を過不足なく表現するために，角に対する鏡映操作を定義し，平面と同様の手順
で壁面寄与成分を計算する．角は 2 つのポリゴンが交わる線分として定義される 2 次元角と，3 つ以上のポリゴンが接する点として定義される 3 次元角に分類される．角と平面に対する鏡映操作は鏡映行列の計算に用b)Refrection

Fig.2 Inconsistency of wall domain around corners occurredin the IERP model

Fig.3 A beam element discretizing a square column member

Fﬂuid
ie2 =

いる単位法線ベクトル n の定義のみが異なる．角 k のni,k は以下のように定義する．

ni,k =

xi − xnear
i,k|xi − xnear
i,k |

(22)

xnear
i,k は角 k 上の粒子 i に対する最近傍点位置であり，3次元角の場合は角 k の位置と一致する．式 (22) の定義
では，壁面境界条件における法線ベクトルが粒子の位置に応じて滑らかかつ連続に変化しており，これは角における壁面境界条件の不連続性が解消されるようなモデリングになっていることを意味する．
c) 非整合領域の体積補正
Fig.2 に示すように，壁領域の形状によっては本来の壁領域と IERP モデルで表現される壁領域に差異が生じる．この誤差を解消するため，IERP モデルでは角に対して計算される壁面寄与成分に次式の体積補正係数 εkを乗じる．
θdεi,k = k(23)θkθk は角 k の流体領域側の角度，θkd は角 k に対する鏡映操作で本来表現されるべき壁領域の角度である．角が2 次元角の場合 θ は平面角，3 次元角の場合 θ は立体角として計算される．詳細な計算方法は文献 [8] を参照されたい．
(2) はり要素に作用する節点流体力の計算
IERP モデルによって計算された壁面に作用する流体
力を，はり要素の等価節点外力に変換する．例として，
Fig.3 のような要素 ie で角柱部材を離散化している場合を考える．要素 ie が表現する部材を構成するポリゴン壁
および角の集合を Wie と定義する．流体粒子 i がポリゴンもしくは角 k ∈ Wie に与える力を fi,k とすると， fi,k の作用点は xnear
i,k である．LBTE の形状関数 N1 (ξ), N2 (ξ) を用いて fi,k を等価節点力に変換すると，要素 ie のローカ
ﬂuid
ル節点 1, 2 に作用する流体力の合力 Fﬂuid
ie1 , Fie2 は，それぞれ次式のように得られる．
∑ ∑{
}
Fﬂuid
=
N1 (ξi,k ) fi,k(24)
ie1k∈Wie i∈N∑ ∑{N2 (ξi,k ) fi,k

}

(25)

k∈Wie i∈Nξi,k は xnear
i,k に対応するはり要素上の正規化座標値である．
また，流体粒子がはり要素 ie に与える流体力に伴って生じるトルクは，
∑ ∑[
{
(
) }]
N1 (ξi,k )ri,k × fi,k − fi,k · eie eie(26)
Tiﬂuid=e1k∈Wie i∈N

Tiﬂuid
=
e2

∑ ∑[

{
(
) }]
N2 (ξi,k )ri,k × fi,k − fi,k · eie eie(27)k∈Wie i∈Nとなる．ただし，ri,k はトルクの腕に相当するベクトル，eie は要素 ie の軸方向単位ベクトルである．

(3) ポリゴン頂点の位置・速度・加速度の更新
次に，ポリゴン頂点の更新計算について述べる．要素 ie のローカル節点 1, 2 それぞれに対応するポリゴン
頂点の集合を Vie1 , Vie2 と定義する．第 n ステップにおける頂点 l の位置は，()
n−1xnl∈Vi = xnie1 + Snie xn−1(28)l∈Vi − xie1
e1

e1

(

n−1
xnl∈Vi = xnie2 + Snie xn−1l∈Vi − xie2

)

e2

e2

(29)

となる．ただし，Snie は n ステップ目における要素 ie の回転行列である．
頂点の速度は，節点速度に軸回りの回転速度成分を足し合わせることで得られる．要素 ie の軸回りの角速度を ωie とすると，頂点 l の速度は，()
vnl∈Vi = vnie1 + ωnie xnl∈Vi − xnie1 × eie
(30)
e1

e1

(

)
vnl∈Vi = vnie2 + ωnie xnl∈Vi − xnie2 × eie
e2

e2

(31)

と計算できる．加速度についても同様である．

4. 開発手法の計算フロー
3 節で記した連成モデルを用い，Conventional Serial
Staggered (CSS) 法 [18] に基づいて ISPH 法と ASI-Gauss
法の分離型弱連成解析を行う．1 つのステップ内における計算の流れを以下に記す．

Table1 Numerical conditions of the free fall analysis

Item
Particle spacing
Eﬀective radius
Slip condition
Time increment
Element division number
Damping coeﬃcients

Fig.4 Schematic view of the free fall problem

1. 式 (28–31) などに従って，ポリゴン頂点の位置・速度・加速度を更新する．
2. 式 (1) を安定化 ISPH 法（式 (4, 6–9)）で解き，流体圧力・流速・壁面に作用する流体力を計算する．
3. 式 (24–27) を用い，はり要素に作用する等価節点
外力を計算する．
4. 式 (10) を ASI-Gauss 法で解き，構造物の変位・速度・加速度を計算する．
5. 数値計算例(1) 部材落下問題
直 方 体 の 容 器 に 溜 め た 水 に 向 かって 角 柱 部 材 を
自 由 落 下 さ せ る 問 題 を 解 く．系 の 概 観 を Fig.4 に
示 す．部 材 は 0.04 m 角 の 正 方 形 断 面 で ，密 度 は200, 400, 600, 800 kg/m3 である．解析条件は Table1 に示すとおりである．
Fig.5 に解析終了時（t = 3.0 s）の可視化結果を示す．
図中の R は流体（水）と部材の密度比で，コンターは圧力を表している．図から，平衡状態において部材が浮く位置が異なっていることが確認できる．ここで，平衡状態における浮力と重力のつり合いから，部材全体の体積に対する水に沈んでいる部分の体積比 V sub /V total は
密度比 R と一致する．解析で得られた体積比 V sub /V totalの時刻歴（Fig.6）から，確かに浮力と重力のつり合いが正しく再現されていることがわかる．
Fig.5 から，平衡状態における部材の向きが密度比 R によって異なることがわかる．具体的には，R = 0.2, 0.8 のケースでは部材の面が下を向く形で静止しており，R =0.4, 0.6 のケースでは部材の辺が下を向いている．この
傾向は数学的・数値的に確認されている結果 [19][20] と整合している．
(2) 弾性ゲート付きダムブレイク問題
紙面の都合上，本稿では割愛し，口頭発表にて内容を説明する．
6. 結言
本研究では，津波-建物連成問題を高精度・高効率に解く手法として，安定化 ISPH 法による自由表面流れ解

Value
0.0050 mm
0.0155 mmfree slip
0.0005 s1
α = 0.0, β = 0.0

析と ASI-Gauss 法による骨組構造解析のカップリングに基づく FSI 解析手法を開発した．粒子とはり要素の連成モデルには，IERP 壁境界モデルをベースとして，はり要素の自由度縮約を考慮した節点外力計算モデル
およびポリゴン更新モデルを適用した．精度検証として角柱落下問題を解いた結果，浮力と重力のつり合いおよび平衡状態における部材の向きが妥当に再現され，十分な精度で連成解析を実現できることが確認された．
今後は，より詳細な精度検証の実施と，実問題への適用に取り組む予定である．
謝辞: 本研究は，JST 次世代研究者挑戦的研究プログ
ラム JPMJSP2124 の支援を受けたものである．ここに記して謝意を表する．
参考文献[1] 国土交通省, 東日本大震災における津波による建築物被害を踏まえた津波避難ビル等の構造上の要件に係る暫定指針, 2011.
[2] Antoci, C. et al., Numerical simulation of ﬂuidstructure interaction by SPH, Coumputer & Structures,Vol. 85, pp. 879–890, 2007.
[3] Yang, Q. et al., Free-surface ﬂow interactions with deformable structures using an SPH-FEM model, Ocean
Eng., Vol. 55, pp. 136–147, 2012.
[4] Asai, M. et al., A stabilized incompressible SPH
method by relaxing the density invariance condition,
J. Applied Mathematics, Vol. 2012, 139583, 2012.
[5] Isobe, D., Progressive Collapse Analysis of Structures: Numerical Codes and Applications, eBook
ISBN: 9780128130421, Elsevier, 2017.
[6] 大村浩之ほか, 壁領域の角を考慮したポリゴン壁境
界モデルの開発およびその ISPH 法への適用, 日本
計算工学会論文集, 2021 巻, p. 20210011, 2021.
[7] Mitsume, N. et al., Explicitly represented polygon
wall boundary model for explicit MPS method, Computational Particle Mechanics, Vol. 2, pp. 73–89,
2015.
[8] Koshizuka, S. and Oka, Y., Moving-particle semiimplicit method for fragmentation of incompressible
ﬂuid, Nuclear Science and Eng., Vol. 123, pp. 421–
434, 1996.
[9] 井元佑介ほか, 安定化 ISPH 法の理論的解釈 II–誤差
評価に基づく安定化係数の最適化–, 土木学会論文集
A2, Vol. 75, pp. 187–194, 2019.

R = 0.2

R = 0.4

R = 0.6

R = 0.8

Fig.5 Visualization of results of the free fall analysis at 3.0 s

Analysis1

R = 0.8

0.8Volume ratio

TheoryR = 0.6

0.6

R = 0.4

0.4

R = 0.2

0.2
0
0

0.5

1

1.5

2

2.5

3

Time (s)Fig.6 Time histories of V sub /V total[10] Monaghan, J. J., Smoothed particle hydrodynamics, J. Computational Physics, Vol. 110, pp. 543–574,
1992.
[11] Tanaka, M. and Matsunaga, T., Stabilization andsmoothing of pressure in MPS method by QuasiCompressibility, J. Computational Physics, Vol. 229,
pp. 4279–4290, 2010.
[12] Tsuruta, N. et al., A short note on Dynamic Stabilization of Moving Particle Semi-impicit method, Computers and Fluids, Vol. 82, pp. 158–164, 2013.
[13] Fung, T. C., Numerical dissipation in time-step integration algorithms for structural dynamic analysis,
Prigress in Structural Eng. and Materials, Vol. 5, pp.
167–180, 2003.
[14] Hughes, T. J. R. et al., A simple and eﬃcient ﬁnite
element for plane bending, Int. J. Numerical Methods
in Eng., Vol. 11, pp. 1529–1543, 1977.
[15] Toi, Y., Shifted Integration technique in onedimensional plastic collapse analysis using linear andcubic ﬁnite element, Int. J. Numerical Methods inEng., Vol. 31, pp. 1537–1552, 1991.

[16] Kawai, T., New element models in discrete structural
analysis, J. Society of Naval Architects of Japan, Vol.
141, pp. 174–180, 1977.
[17] Toi, Y. and Isobe, D., Adaptively shifted integration technique for ﬁnite element collapse analysis offramed structures, Int. J. Numerical Methods in Eng.,
Vol. 36, pp. 2323–2339, 1993.
[18] Farhat, C. and Lesoinne, M., Two eﬃcient staggered
algorithms for the serial and parallel solution of threedimensional nonlinear transient aeroelastic problems,Computer Methods in Applied Mechanics and Eng.,
Vol. 182, pp. 499–515, 2000.
[19] Gilbert, E. N., How things ﬂoat, The American Math.
Monthly, Vol. 98, pp. 201–206, 1991.
[20] Morikawa, D. S. and Asai, M., Coupling total Lagrangian SPH-EISPH for ﬂuid-structure interaction
with large deformed hyperelastic solid bodies, Computer Methods in Applied Mechanics and Eng., Vol.
381, 113832, 2021.