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

日本計算工学会

PMS を用いた SPH-梁要素カップリングに基づく流体構造連成解析手法の開発
Development of a FSI analysis method based oncoupled SPH-beam element using PMS
大村浩之 1) 辻勲平 2) 浅井光輝 3) 磯部大吾郎 4)
Hiroyuki Omura, Kumpei Tsuji, Mitsuteru Asai and Daigoro Isobe
1)

博 (工) 防災科学技術研究所 特別研究員 (E-mail: homura@bosai.go.jp)
2)

博 (工) 東北大学 助教 (E-mail: kumpei.tsuji.e1@tohoku.ac.jp)

3)

博 (工) 九州大学 准教授 (E-mail: asai@doc.kyushu-u.ac.jp)
4)

博 (工) 筑波大学 教授 (E-mail: isobe@kz.tsukuba.ac.jp)

For interaction problems between free surface ﬂow and frame structures, we develop a separated FSI analysis
method based on the coupling of SPH and beam element FEM. Speciﬁcally, we propose that representation
scheme of member shapes with particles corresponding to beam elements, ﬂuid structure interaction method
using PMS (Passively Moving Solid), and modiﬁcation of particle volume when a beam element undergoes
large deformation. The validity of the proposed method is shown by solving static and dynamic ﬂuidstructure interaction problems.
Key Words : FSI analysis, SPH, FEM, Beam element, PMS

1. 緒言
本研究では，津波によって建物が流失する流体構造連成（FSI）問題のシミュレーション実現に向けて，SPHと梁要素の有限要素法（FEM）をカップリングさせた分離型 FSI 解析手法を構築する．粒子法には非圧縮性流れのための Incompressible SPH（ISPH）法 [1] を採用し，FEM には骨組構造の崩壊挙動などの強非線形問題を高精度かつ高効率に解くことができる Adaptively ShiftedIntegration（ASI）-Gauss 法 [2] を用いる．これにより，自由表面を伴う流れを受ける建物が倒壊し移流するプロセスを再現する．
粒子法と FEM のカップリングに関する研究例（例えば [3]）は数多く存在するものの，そのほとんどは流体から構造への連成計算において，構造物の境界を壁粒子もしくはポリゴンなどで表現し，境界面に作用する流体力を表面積分するように定式化されている．その場合，連成面の力のつり合いを満たすためには，安定化項を含む 0 次精度の圧力勾配モデルしか適用できないという問題点がある．一方で，境界面上の圧力を表面積分して流体力を得ることも可能だが，この方法では連成面上の力のつり合いが保証されない．さらに，上記の方法では構造物表面に壁面 Nuimann 境界条件を課す必要があるため，構造物が壁面に近接する場合など，1 つの流体粒子に対して異なる方向の Nuimann 条件が同時に課されると圧力計算が破綻する．
以上の課題を解決するために，本研究では Koshizuka
et al., 1998[4] が提案した Passively Massive Solid（PMS）モデルを用いて連成計算を行う．PMS は固体領域も流体と同様に粒子で離散化し，圧力および流体力は固体粒子も仮想的な流体として扱い流体粒子と同様に計算し，固体粒子の速度更新には構造解析で得られた速度を用いる方法である．また，構造物に作用する外力は各固体粒子に作用する流体力の体積積分で得られる．そのため，連成面の力のつり合いは保証され，任意の空間微分モデルを適用できる．かつ，連成面の Nuimann
境界条件は流体粒子と固体粒子の間に課される非圧縮条件で近似されるため，境界条件に起因する計算の破綻も起きない．
これまで，PMS モデルは粒子法で流体剛体連成問題を解くための手法として開発・適用されている [5] が，本研究では流体と骨組構造の連成解析への適用を試みる．そのために，梁要素の部材形状を粒子で離散表現し，要素の変形に応じて固体粒子の体積を補正する方法を提案する．

2. SPH 法に基づく非圧縮性自由表面流れ解析(1) 基礎式と空間微分モデル
非圧縮性流れの Navier-Stokes 方程式は，
Dv
1= − ∇p + ν∇2 v + g
Dtρ

(1)

である．ただし，t, v, ρ, p, g はそれぞれ時間，速度，密度，圧力，重力加速度である．これを SPH および射影法に基づき離散化すると，粒子 i に関して以下の時間発展式が導かれる．
(⟨⟩n)
vi∗ = vni + ν∇2 v + g ∆t(2)i

1vn+1= v∗i − ⟨∇p⟩n+1∆tiiρi⟨⟩n+1 ρi
ρi − ⟨ρ⟩ni∇2 p=⟨∇ · v⟩∗i + αstabi∆t∆t2(3)(4)

ここで，n は任意の時間ステップ，∗ は n と n + 1 ステップの中間における仮想ステップを表す．また，⟨⟩ は粒子で離散化された量であることを表す．式 (4) の右辺第二項は体積保存性を付与するためのソース成分であり，αstab は 1 より十分小さい正の実数である [1]．
本稿では，勾配には Bonet and Lok, 1999[6] の 1 次精度モデル，Laplacian には二次精度モデルとして Asai etal., 2023[7] の full inverse モデルをそれぞれ用いる．
勾配モデルは，∑V j ϕi j ∇wi j(5)⟨∇ϕ⟩i = Lij∈Pi−1∑
Li =  V j ∇wi j ⊗ xi j (6)j∈Piである．ここで，Pi は粒子 i の影響半径内に存在する粒子集合，wi j は粒子 i, j 間のカーネル関数である．また，
ϕi j = ϕ j − ϕi である．発散 ⟨∇ · ϕ⟩i も上式と同様に計算
˜ i j ≡ Li ∇wi j と略記する．
できる．以下，∇w
Laplacian モデルは，次式を解くことで得られる．
 ⟨ ∂2 ϕ ⟩  ⟨ (∂x1 )2 ⟩i  ∂2 ϕ ⟨ (∂x2 )2 ⟩i  ∂2 ϕ (∂x )2 ∑{) }( ⟨ 3 i⟩  = 2G−1V jCi j ϕi j − xi j · ⟨∇ϕ⟩i si j (7)


i
∂2 ϕ
2 ∂x1 ∂x2 j∈Pi ⟨⟩i 
∂2 ϕ2 ∂x2 ∂x3  ⟨ 2 ⟩i 
2 ∂ϕ 
∂x3 ∂x1iCi j =Gi =∑˜ ijxi j · ∇w
∥xi j ∥4(8)V j C i j t i j ⊗ si j

(9)

j∈Pi

 2 
 xi j,1  2  xi j,2  2  xsi j =  i j,3  xi j,1 xi j,2  xi j,2 xi j,3 xi j,3 xi j,1∑
2 ˜xi2j,1 − xi j · k∈Pi Vk xik,1∇wik


∑


22 ˜xi j,2 − xi j · k∈Pi Vk xik,2 ∇wik

∑


2 ˜xi2j,3 − xi j · k∈Pi Vk xik,3∇wik
t i j = ∑˜ ik 
 xi j,1 xi j,2 − xi j · k∈Pi Vk xik,1 xik,2 ∇w∑
˜ ik 
 xi j,2 xi j,3 − xi j · k∈Pi Vk xik,2 xik,3 ∇w
∑
˜ ik xi j,3 xi j,1 − xi j · k∈Pi Vk xik,3 xik,1 ∇w

IERP モデルでは，壁面を多角形ポリゴンで表現し，壁面ポリゴン k に関して流体粒子 i を鏡映したものを鏡映粒子 ik とする．鏡映粒子 ik が受ける流体からの寄与をさらに鏡映することで，流体粒子 i が壁 k から受ける寄与を得ることができる．例えば勾配の壁寄与 ⟨∇ϕ⟩Wiは， ∑∑ ()

WkkRi,k V
ϕ−ϕ∇w⟨∇ϕ⟩i = Li jiji j (12)j∈Pikk∈Wiと書ける．Wi は粒子 i の影響半径内に存在する壁ポリゴンの集合，Ri,k は粒子 i をポリゴン k に関して鏡映する行列である．ここでは割愛するが，発散や Laplacianの壁寄与も同様に導かれる．
鏡映粒子が持つ速度および圧力は，壁面上における速度の Dirichlet 条件および圧力の Nuimann 条件に従い，流体粒子の物理用をマッピングすることで得る．すなわち， Ri,k vi (free-slip)
vik = (13)−vi(no-slip)

pik = pi − ρg · xkiとする． xki はポリゴン k から粒子 i に向かう最短経路ベクトルである．

(3) 粒子再配置処理
一般に，SPH 法では粒子の空間的な配列が乱れるほど微分モデルの計算精度が低下するため，粒子の配列を逐次整える処理（粒子再配置処理）が必要になる．ここでは Optimized PS（OPS）[9][10] に基づく粒子再配置処理を導入する．
3. ASI-Gauss 法に基づく骨組構造解析(1) 基礎式と線形 Timoshenko 梁要素
Newton-Raphson 法に基づき導出された動的平衡方程
式は，
Ku + Cu̇ + Mü = f − (q + Cv + Ma)

(10)

(11)

(2) 壁境界モデル
本稿では，Improved Explicitly Represented Polygon（IERP）モデル [8] を用いて壁面境界条件の処理を行う．ここでの壁面境界は流体構造連成境界とは区別しており，空間上に固定された壁面を指す．

(14)(15)と書ける．K, C, M, f , q はそれぞれ接線剛性行列，減衰行列，質量行列，外力および内力ベクトルである．u, v, aはそれぞれ変位，速度，加速度ベクトルである．式 (15)に HHT 法 [11] を適用すると，n + 1 ステップの i イタ
レーションにおける解くべき方程式として以下が得られる．
)(γ n+1,i1n+1,in+1,iM∆uK+C+β∆tβ∆t2(16)()
n+1n+1,i n+1,in+1,i n+1,in+1,i= f− Cv+Ma+qここで，∆u は変位修正量である．また，β, γ は積分パラメータであり，HHT 法に基づき

β=1(1 − α)24(17)1−α2(18)γ=と表せる．α は −1/3 ≤ α ≤ 0 を満たす数値減衰パラメータであり，値が小さいほど減衰効果が強くなる．
n + 1 ステップにおける変位，速度，加速度の初期値は次式の通り求められる．

un+1,0 = 0()()γ nγ nn+1,0v= 1−v + ∆t 1 −aβ2β()1 n1 nn+1,0a=−v + 1−aβ∆t2βSection
𝑆(𝜉)

Γ𝒙(𝜉)

(19)
−1𝜉(20)1Particle representation(21)

𝕊𝜉𝔾𝒙(𝜉)

速度修正量 ∆v，加速度修正量 ∆a は，∆v =γ∆uβ∆t1∆a =∆uβ∆t2(22)𝜉1(23)である．
M, K, f , q の計算には，線形 Timoshenko 梁要素 [12]（Linear Timoshenko Beam Element, LTBE）を用いる．
LTBE は 2 節点 12 自由度の梁要素であり，全ての変位成分を同一の形状関数で補間する．形状関数は以下に示す 1 次 Lagrange 補間関数である．

1N1 (ξ) = (1 − ξ)
2(24)1(1 + ξ)2(25)

N2 (ξ) =−1

ξ は自然要素座標の値であり，定義域は −1 ≤ ξ ≤ 1 である．また，C は Rayleigh 減衰モデルを用いて求める．
(2) ASI-Gauss 法
ASI-Gauss 法 [2] は，LTBE の降伏を塑性ヒンジモデルで表す際に，塑性ヒンジの発生箇所を断面力の外挿計算で予測したうえで，その位置の断面力を直接評価するように数値積分点をシフトするスキームである．これにより，最小限の要素分割で塑性，破壊を含む強非線形性の挙動を高精度に解くことができる．さらに，LTBEを 2 要素 1 サブセットで扱い，サブセット全体で 2 点積分を行うように各要素の積分点を配置することで，弾性域における曲げ変形解の精度を向上させている．
(3) 粒子ベース接触モデル
本研究では，部材表面を流体と構造の連成境界面として扱うため，部材形状を考慮しない一般的な梁要素の接触計算手法を用いると，接触面近傍の連成計算が破綻する．この問題を解消するため，本研究では著者が開発した粒子ベース接触モデル [13] を用い，部材形状を補間する粒子に接触拘束条件を課すことで部材表面の接触現象を精緻に解く．
4. PMS を用いた FSI 解析(1) 梁部材形状の粒子補間
等断面梁要素の立体形状を，断面を表す有限平面領
域 S (ξ) を部材軸方向に補間することで表すことを考え

図–1 Representation of a beam element with particlesる．Timoshenko 梁理論における断面剛の仮定から，時刻 t から t + ∆t にかけて，断面内の任意の位置 x ∈ S (ξ)
は次式のように更新される．

xt+∆t = Pt+∆t (ξ)xt 1
P(ξ) ≈  ∆θz (ξ)
−∆θy (ξ)−∆θz (ξ)
1∆θ x (ξ)(26)
∆θy (ξ) −∆θ x (ξ)

1(27)ここで， P は断面の回転行列，∆θ は増分回転変位である．なお，式 (27) は 2 次以上の項を無視して近似的に得ている．
図-1 のように，梁要素の部材形状を有限個の粒子（固体粒子）で離散表現し，その集合を G とする．断面領域 S (ξ) に対応する粒子の部分集合を Sξ と定義する．Gおよび Sξ (−1 ≤ ξ ≤ 1) は時間発展によらず変化しないものとする．
粒子 i ∈ St+∆tの位置 xi は，要素上の ∆ut+∆t (ξ) と式ξ(26) を用いて次式のように得られる．

xt+∆t= xti + ∆ut+∆tii()t+∆tt+∆tt+∆t∆ui = ∆u (ξ) + P (ξ) − I rti

(28)
(29)

ここで， ri は S (ξ) の中心（要素軸上）に対する粒子 iの相対位置ベクトルである．

(2) PMS に基づくカップリング計算フロー
PMS に基づく流体構造カップリングの計算手順は以下に示す通りである．
1. 固体粒子も非圧縮流体として扱い，全粒子に関して式 (2) における粘性項を求める．
2. 流体粒子のみ式 (2) に従い速度を更新する．
3. 固体粒子も非圧縮流体として扱い，全粒子に関して式 (4) に従い圧力場，式 (3) における圧力勾配項を求める．
4. 流体粒子のみ式 (3) に従い速度と位置を更新する．
5. 1. および 3. それぞれで求めた固体粒子に作用する粘性力および圧力勾配力を固体領域全体で積分することで，固体に作用する流体力を得る．

Unit: [m]6. 5. で求めた流体力を外力として構造解析を行い，固体粒子の速度と位置を更新する．
梁要素上の ξ に対応する固体粒子 i が受ける流体力
fiﬂuid は，(⟨)⟩
1fiﬂuid = mi ν∇2 v − ⟨∇p⟩i
(30)
i
ρi

Loading
0.25

0.20と計算される． fiﬂuid を用い，梁要素の局所節点 1, 2 それぞれに作用する流体力 f および流体トルク t は次式のように計算できる．
∑f1 =
N1 (ξ) fiﬂuid
(31)
i∈G

f2 =
t1 =

∑i∈G

t2 =

∑∑

N2 (ξ) fiﬂuid

(32)

{}
N1 (ξ) ri × (I − e ⊗ e) fiﬂuid(33){}N2 (ξ) ri × (I − e ⊗ e) fiﬂuid(34)i∈Gi∈GX0.40

図–2 Initial conﬁguration of the hydrostatic problem
表–1 Fluid analysis conditions for the hydrostatic problem

Time increment [s]Particle spacing [m]Stabilization parameter αstab
Kinematic viscosity [m2 /s]
Density [kg/m3 ]5.0 × 10−4
5.0 × 10−3
5.0 × 10−3
1.0 × 10−6
1.0 × 103

表–2 Structural analysis conditions for the hydrostatic prob-lem

ここで，e は要素軸方向の単位ベクトルである．

(3) 梁要素の変形に応じた粒子体積補正
梁要素の変形に伴って固体粒子が変位すると固体粒子の空間分布に粗密が生じる．この状態のまま流体解析を行うと，SPH における微分モデルの計算精度が低下する．また，変形が進行し分布が疎になりすぎると，固体領域内に流体粒子が侵入する可能性もある．そこで，本研究では要素の変形，すなわち体積変化に固体粒子の体積を整合させるために，梁要素の体積変化率に基づく固体粒子の体積補正計算を導入する．
変形勾配テンソル F は次式で定義される．
F= I+∇⊗u(35)

いま，梁要素は粒子で離散表現されているため，式 (35)における微分項は SPH に従って計算できる．すなわち，固体粒子 i 上の増分変形勾配テンソル ⟨∆F⟩i は次式を用いて計算できる．
∑˜ i j ⊗ ∆ui j
V j ∇w(36)⟨∆F⟩i =j∈(Pi ∩G)˜ i j の計算に使用すなお，ここでの修正カーネル勾配 ∇w
る Li は，式 (6) における近傍粒子集合 Pi を Pi ∩ G に置き換えて計算する．
体積変化率は変形勾配テンソルの行列式として定義される．時刻 t + ∆t における固体粒子 i の体積 Vit+∆t は次式のように求められる．

00.40ZVit+∆t = κit+∆t V 0(37)κit+∆t = det ⟨∆F⟩ti κit(38)

ここで，V は初期粒子体積であり，粒子径 d を用いてV 0 = d3 で与えられる．
式 (38) に従って逐次的に体積変化率を更新する場合，粒子で近似計算することに起因する誤差の蓄積が生じTime increment [s]Number of element divisions
Parameter of HHT method αWidth [m]Young’s modulus [N/m2 ]Poisson’s ratio
Density [kg/m3 ]5.0 × 10−44−1/3
5.0 × 10−2
7.93 × 107
0.400.491 × 103

る．実際，梁要素が剛体運動する場合でも，時間経過に伴って体積変化率が増大することを予備計算によって確認している．ここでは，梁要素を構成する固体粒子集合 G における平均体積変化率はその要素の軸方向ひずみ εzz に対応すると仮定することで，誤差蓄積を回避するための計算を導入する．すなわち，体積変化率 κiを次式の通り修正し，それを体積補正に用いる．

κi B|G| (1 + εzz )∑κii∈G κi(39)

5. 計算例(1) 曲げ変形を受ける梁部材まわりの静水圧はじめに，プログラムの基本的な動作ならびに体積補正計算の効果を検証するために，曲げ変形を受ける梁部材が設置された容器内の静水圧問題を解く．対象とする系の概観を図-2 に示す．部材の断面形状は正方形で，容器の中央に片端固定で設置する．
部材の曲げ変形は先端に集中横荷重を与えることで生じさせる．流体解析と構造解析それぞれの解析条件を表-1 および表-2 に示す．滑り条件は free-slip とし，構造物の物理減衰は導入していない．体積変化率を常に1 とするケースと，式 (37) による体積補正を行うケースの 2 通りで計算を実施する．
解析で得られた可視化結果，流速分布，圧力分布をそれぞれ図-3，図-4，図-5 に示す．図-3 は Y = 0.2 [m]の断面を表しており，流体のコンターは速さ，構造のコンターは体積変化率である．図-4，図-5 の横軸は水Without volume modification

ZYX
With volume modification

図–3 Visualization result of the hydrostatic problem

図–5 Pressure distribution obtained from simulation of thehydrostatic problem
Unit: [m]Impact1.91.5

図–4 Velocity distribution obtained from simulation of thehydrostatic problem

深である．可視化結果と流速分布から，固体粒子の体積を補正していない場合，部材の固定端，すなわち変形角が最も大きい位置で流速に非物理的な乱れが生じていることが分かる．一方で体積を補正した場合では，部材の周辺で若干のノイズは認められるものの，上述のケースに対して大きく乱れが解消されていることが確認できる．圧力分布についても，流速ほど顕著ではないが，体積補正によって攪乱が低減されている．

(2) 水中における片持ち梁の自由振動
流体を伴う構造物の動的挙動に対する再現性を検証するために，静止させた水の中における片持ち梁の自由振動問題を解く．比較対象として，Liu et al., 2014[14]の実験結果を参考とする．対象とする系の概観を図-6に示す．部材断面の形状は正方形で，片端固定で容器内中央に部材を設置している．実験では容器の寸法は図-6 のものよりさらに大きいが，ここでは計算時間の都合から縮小している．梁部材の先端に衝撃荷重を与えることで自由振動を起こす．
流体解析と構造解析それぞれの解析条件を表-3 およZ

0.09
X0.3

図–6 Initial conﬁguration of the free-vibration problem表–3 Fluid analysis conditions for the free-vibration problem

Time increment [s]Particle spacing [m]Stabilization parameter αstab
Kinematic viscosity [m2 /s]
Density [kg/m3 ]5.0 × 10−4
3.0 × 10−3
5.0 × 10−3
0.903 × 10−6
0.997 × 103

び表-4 に示す．側壁の影響を極力小さくするように，滑り条件は free-slip としている．実験で測定された空気中での振動特性と一致するように，部材の密度と減衰定数をキャリブレーションしている．
解析で得られた部材の変位が最大となる時刻 0.14 [s]における可視化結果を図-7，空気中・水中それぞれにおける 1 次の固有周波数と減衰比を表-5 に示す．まず可あることを示した．
表–4 Structural analysis conditions for the free-vibration

problemTime increment [s]Number of element divisions
Parameter of HHT method αWidth [m]Young’s modulus [N/m2 ]Poisson’s ratio
Density [kg/m3 ]5.0 × 10−4
12−1/55.0 × 10−22.30 × 1011
0.25
7.80 × 103

ZYX

図–7 Visualization result of the free-vibration problem
表–5 The ﬁrst frequencies and damping ratios obtained fromexperiment and simulation

Frequency in air [Hz]Damping ratio in air
Frequency in water [Hz]Damping ratio in water

Experiment[14]7.260.00187.070.0075

Simulation
7.270.0018
7.040.0101

視化結果から，流速は部材の周辺にのみ分布しており，側壁の影響は極めて小さいことが確認できる．水中における固有周波数と減衰比は，実験結果に対してやや過大に評価されているが，これは固体と流体間の滑りを正確に再現できていないためと思われる．PMS における滑り条件の付与は今後の課題としたい．

(3) 段波と漂流物の衝突を受ける骨組構造物の流失提案手法の有用性を示すために，骨組構造物が段波および漂流物の衝突を受けることで流失する過程のシミュレーションを実施する．結果は口頭発表にて詳述する．
6. 結言
本研究では，PMS を用いた SPH-梁要素 FEM カップリングに基づく流体構造連成解析手法を開発した．PMSを梁要素とのカップリングに適用するために，梁要素の立体形状を考慮した粒子離散表現を定式化した．また，
PMS を流体構造連成に適用するために，構造物の変形に伴う粒子体積補正計算を提案した．提案手法を用い，曲げ変形を受ける梁部材まわりの静水圧問題，ならびに水中における片持ち梁の自由振動問題を解き，理論値もしくは実験値と整合する結果が得られた．また，段波と漂流物の衝突を受ける骨組構造物の流失問題を解くことで，提案手法が津波被害予測技術として有用で

参考文献[1] M. Asai et al., A Stabilized Incompressible SPH
Method by Relaxing the Density Invariance Condition.
J. Appl. Math. 2012 (2012).
[2] K. M. Lynn and D. Isobe, Finite element code for impact collapse problems of framed structures. Int. J. Numer. Methods Eng. 69, 2538–2563 (2007).
[3] S. W. Attaway et al., Coupling of smooth particle hydrodynamics with the ﬁnite element method. Nucl.
Eng. Des. 150, 199–205 (1994).
[4] S. Koshizuka et al., Numerical analysis of breaking
waves using the moving particle semi-implicit method.
Int. J. Forecast. 26, 751–769 (1998).
[5] H. Negishi et al., Fluid-rigid body coupling simulations with the passively moving solid model based ona physically consistent particle method. Phys. Fluids.
36 (2024).
[6] J. Bonet and T. S. L. Lok, Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations. Comput. Methods Appl. Mech.
Eng. 180, 97–115 (1999).
[7] M. Asai et al., A class of second-derivatives in theSmoothed Particle Hydrodynamics with 2nd-order accuracy and its application to incompressible ﬂow simulations. Comput. Methods Appl. Mech. Eng. 415,116203 (2023).
[8] 大村 浩之 ほか, 壁領域の角を考慮したポリゴン壁境界モデルの開発およびその ISPH 法への適用. 日本計算工学会論文集 2021, 20210011 (2021).
[9] S. J. Lind et al., A generalised diﬀusion-based algorithm for stability and validations for impulsive ﬂows
and propagating waves. J. Comput. Phys. 231, 1499–1523 (2012).
[10] A. Khayyer et al., Comparative study on accuracy
and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting scheme in ISPH context. J. Comput. Phys.
332, 236–256 (2017).
[11] H. M. Hilber et al., Improved numerical dissipation
for time integration algorithms in structural dynamics.
Earthq. Eng. Struct. Dyn. 5, 283–292 (1977).
[12] T. J. R. Hughes et al., A simple and eﬃcient ﬁnite
element for plate bending. Int. J. Numer. Methods Eng.
11, 1529–1543 (1977).
[13] 大村 浩之, 梁要素を用いた FEM のための粒子
ベース接触計算手法. 日本計算工学会論文集 2024,20241004 (2024).
[14] F. Liu et al., Added mass matrix estimation of beams
partially immersed in water using measured dynamicresponses. J. Sound Vib. 333, 5004–5017 (2014).