A-04-03

第26回計算工学講演会

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

計算工学会

角付近の壁面寄与を考慮した改良 ERP 壁境界モデルの開発とその ISPH 法への適用
Development of Improved Explicitly Represented Polygon Modelwith Consideration of Contribution of Wall near Edges
and Application to ISPH 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 paper, we develop polygon wall boundary model with consideration of edges of wall domain thatis suitable for semi-implicit particle method, such as Incompressible Smoothed Particle Hydrodynamics(ISPH) method. The polygon model is based on Explicitly Represented Polygon (ERP) model, and deﬁnitionof reﬂecting operation for surface polygon is extended in order to eliminate error occurred near edges in the
ERP model. Besides, volume of wall domain expressed near edges are modiﬁed using an approximate
polynomial based on planer angle and solid angle of edges to improve accuracy of the calculation. We
introduce the improved ERP model to ISPH method, and solve a three-dimensional hydrostatic problem as anumerical example. According to the results, it is conﬁrmed that the improved ERP model can obtain moreaccurate and smoothed distributions of pressure than the previous ERP model.
Key Words : Particle Method, Incompressible Smoothed Particle Hydrodynamics Method, Polygon
Wall Boundary Model, Solid Angle

1. 序論
粒子法に基づく流体解析手法は現在までに数多く開発されており [1, 2]，近年では流体構造連成（FluidStructure Interaction, FSI）解析にも利用されている [3,4]．FSI の流体解析スキームとして粒子法を用いる場合，
構造解析側へ受け渡す外力として，壁面に作用する流体力を計算する必要がある．そのためには，壁面境界
が流体領域に及ぼす影響をモデル化する必要がある．
粒子法における壁境界モデルは，空間上に固定された粒子で壁面を表現するモデル（壁粒子モデル）と，壁面境界の位置にポリゴンメッシュを配置し，メッシュに作用する流体力を直接計算するモデル（ポリゴンモ
デル）に大別される．特に，ポリゴンモデルは壁面形状を滑らかに再現できるため，幾何学的な整合性が担保されることから FSI への適用性が高いといえる．その中でも，Mitsume et al. による Explicitly Represented
Polygon 壁境界モデル [5] は，壁面境界上での力学的な物理量が整合的に計算できるため，FSI に適したポリゴンモデルとして知られている．
粒子法には圧力を陽的に求める手法と陰的に求める手法の 2 つがあるが，FSI へ組み込むことを考えた場
合，圧力解が陽解法と比べてより高精度に求まる半陰解法を使用することが望ましい．ところが，ERP 壁境界モデルでは角付近で生じる局所的な計算誤差が生じることが指摘されており [5]，圧力を陰的に計算する場© 一般社団法人 日本計算工学会

合，この局所的な誤差は系全体に悪影響を及ぼすことが予想される．したがって，ERP 壁境界モデルを半陰解法に適用するためには，上記の角付近における計算誤差を解消する必要がある．
以上を踏まえ，本研究では，精緻な FSI 解析実現に向けた基礎的検討として，従来の ERP 壁境界モデルに対して，角付近における壁面寄与の計算精度を改善した
改良 ERP 壁境界モデルを開発する．具体的な方法としては，従来モデルでは平面に対してのみ定義されていた壁面寄与の計算を角に対しても定義する．そのうえ
で，角付近の壁領域評価の精度を保証するために，角に対して表現される壁領域の体積を補正する処理を導入する．この際，計算効率を損なわないように，角の幾
何学的な情報に基づく単純な代数式を用いて近似的に体積補正を行う．また，改良 ERP 壁境界モデルを半陰的粒子法の一種である Incompressible Smoothed Particle
Hydrodynamics (ISPH) 法に適用し，3 次元の静水圧問
題を解くことで，提案モデルの妥当性と有用性を検証する．

2. ISPH 法の概要
(1) 支配方程式
非圧縮性流体の運動は次式の質量と運動量の保存則
によって支配される．

- A-04-03 -

∇·u=0

(1)

A-04-03

第26回計算工学講演会(4) SPH 法に基づく空間離散化スキーム
SPH 法は，空間中の任意の点 x 上で定義される任意のスカラー関数 ϕ(x) が，次式のように周囲の点 ξ 上の
物理量 ϕ(ξ) を参照してカーネル近似できることを基礎としている．
∫
ϕ(x) =
W(|x − ξ|, h)ϕ(ξ)dξ(7)

ParticleΠ

図–1 粒子法に基づく離散化

1Du= − ∇p + ν∇2 u + g
Dt
ρ

(2)

ただし，ρ は密度，u は流速， p は圧力， g は重力加速度を表す．式 (1)(2) それぞれにおける時間微分項を離散化することで，速度と圧力を求めるための方程式が得られる．

(2) 時間積分スキーム
次に，速度と圧力の更新手順を示す．まず，式 (2) 左辺の時間微分項を以下に示す前進差分の形で離散化する．
Du un+1 − u∗ u∗ − un
≈
+
Dt
∆t
∆t

(3)

ただし，u∗ は中間状態における速度である．また，∆tは解析の時間増分である．ここで，中間状態はナビエ
ストークス方程式（式 (2)）の右辺第一項である圧力勾配項を無視した仮の状態であると定義する．すなわち，u∗ は次式を解くことで求まるものとする．

u∗ − un= ν∇2 un + g∆t(4)

なお，u∗ は速度予測子と呼ばれる．
続いて，速度予測子と圧力勾配項から n + 1 ステップにおける速度，すなわち速度修正子を計算する．なお，
圧力勾配項を計算するためには圧力を求めておく必要があるが，その方法は次項で述べる．速度修正子は次式を解くことで求まる．

un+1 − u∗1= − ∇pn+1
∆t
ρ

(5)

ここで，W はカーネル関数であり，特定の条件を満たす重み関数である．本研究ではカーネル関数として 5
次スプライン関数を用いる．h は積分の有効半径で，粒子法においては影響半径と呼ばれる．
図-1 に示すように，物体領域を有限個の粒子（計算点）で離散化し，各粒子の質量 m と密度 ρ を定義すれ
ば，式 (7) は次式で表す重み付き総和の形に近似できる．
∑ mj
ϕ(xi ) ≈ ⟨ϕi ⟩ =W(|ri j |, h)ϕ j
(8)
ρj
j∈P
i

ただし，⟨ ⟩ は SPH 法による近似値を表す．i および j は
粒子番号を示し，ri j = xi − x j は粒子 j に対する粒子 i の相対位置ベクトルである．また，以降は W(|ri j |, h) = Wi jと略記する．Pi は粒子 i の影響半径内に存在する粒子の集合を表す．
オリジナルの SPH 法では，圧力の勾配，速度・圧力のラプラシアン，速度の発散の離散化は式 (8) に空間微分を適用することでそれぞれ導かれる [1]．しかし，その定式化では計算精度や数値安定性に問題があるため，
現在では様々な改良式が提案されている．本研究では，
Asai et al.[6] にならい，以下に示す離散化式を採用する．
• 圧力の勾配
 



∑


  p j pi m j  2 + 2  ∇Wi j (9)⟨∇pi ⟩ = ρi





ρ
ρ
j∈Pi

i

• 速度のラプラシアン)}⟨⟩ 2 ∑ { ( ri j · ∇Wi j ) (∇2 ui =
mj
u
−
u
i
j
ρi j∈P|ri j |2 + η2

(10)

i

• 圧力のラプラシアン⟨⟩ 2 ∑ { ( ri j · ∇Wi j ) ()}
2∇ pi =
mjpi − p j
ρi j∈P|ri j |2 + η2

(11)

i

• 速度の発散⟨∇ · ui ⟩ =(3) 圧力ポアソン方程式
圧力解を求めるためのポアソン方程式を次式に示す．

j})1 ∑{ (
m j u j − ui · ∇Wi j
ρi j∈P

(12)

i

(6)

ただし，式 (10)(11) の η はゼロ除算を防ぐための安定化パラメータであり，文献 [6] を参考に η = 0.01h と設定している．

右辺のソース項には Asai et al.[6] の非圧縮条件緩和型
ISPH 法で提案されているものを使用する．α は密度一定条件，すなわち非圧縮条件を緩和するための係数であり，時空間の解像度に応じて値が設定される．

3. 改良 ERP 壁境界モデル(1) 鏡映行列と鏡映粒子の定義
図-2 に示すように，流体粒子 i がポリゴン k で表現される壁領域から受ける寄与 ϕwall
i,k は，近傍流体粒子群∇2 pn+1 ≈ −© 一般社団法人 日本計算工学会

ρρ −ρ
∇ · u∗ + α
∆t
∆t2
0

0

n

- A-04-03 -A-04-03

第26回計算工学講演会

Fluid domain
Wall domain

Polygon

図–2 流体粒子 i がポリゴン k で表現される壁領域から受ける

図–5 角周りで生じる非整合領域（赤色の網掛け範囲）と鏡映

寄与 ϕwalli,k

操作によって表現される壁領域（灰色の粒子群）

計算される．

ni,k =

xi − xnear
i,k
=
near
|x
−
x
|rwall
|
i
i,k |
i,k
rwall
i,k

xnear
i,k は粒子 i に対するポリゴン k 上の最近傍点の位置ベクトルである．なお，rwall
i,k はポリゴン k から粒子 i に向かう法線ベクトルとしている．つまり，あるポリゴ
ン k は， rwall
i,k が面と直交する領域に存在する粒子にのみ寄与するとしている．また，ポリゴン k に対して生成される鏡映粒子 i′ の位置ベクトル xi′ ,k は次式のように計算できる．

図–3 鏡映粒子 i′ が流体領域から受ける寄与 ϕﬂuidi′ ,kxi′ ,k = xi − 2rwalli,k

3-D edge
2-D edge

図–4 2 次元角と 3 次元角j ∈ Pi をポリゴン k に関して鏡映変換した粒子群 j′ ∈ Piとの相互作用で計算できる．
′
このとき，ϕwall
i,k は，鏡映粒子 i がその近傍の流体粒
ﬂuid
子群 j ∈ Pi′ から受ける寄与 ϕi′ ,k （図-3）と鏡映対称の関係になっている．すなわち，粒子 i に対するポリゴ
ン k に関する鏡映変換の表現行列（以下，鏡映行列）を
Ri,k とすれば，ϕwall
i,k は次式のように表せる．
ﬂuid
ϕwall
i,k = Ri,k ϕi′ ,k(13)

ここで，鏡映行列 Ri,k は以下のように定義される．

Ri,k = E − ni,k ⊗ ni,k

(14)

ただし，E は単位行列，ni,k はポリゴン k の面上で定義される外向き単位法線ベクトルであり，次式のように© 一般社団法人 日本計算工学会

(15)

(16)

(2) 角に対する鏡映操作
前項で示したように，ERP 壁境界モデルでは流体領域をポリゴン面に関して鏡映することで壁領域を表現する（図-2, 3 参照）．したがって，角付近ではその近傍のポリゴン面それぞれに対して鏡映操作を行うため，角の周りにおいて表現している壁領域に欠損もしくは重複が生じる（以下，これを非整合領域と呼称する）．
改良 ERP 壁境界モデルでは，この非整合領域を解消することで，角付近の壁面寄与の計算精度を改善する．
はじめに，角を図-4 に示すような「2 次元角」と「3次元角」に分類する．2 次元角は隙間なく隣接する 2 つポリゴン面が接する線分，3 次元角は隙間なく隣接する
3 つ以上のポリゴン面が接する点として定義する．
a) 2 次元角に対する鏡映操作
2 次元角の番号を k∗ として， rwalli,k∗ は通常のポリゴン面に対する場合（式 (15) 参照）と同様に計算する．ただ∗し，xneari,k∗ は 2 次元角 k を形成する線分上における粒子i との最近傍点とする． rwalli,k∗ が決定すれば，式 (13)(16)を用いて角 k∗ に対する鏡映行列 Ri,k∗ と鏡映粒子の位置
xi′ ,k∗ を計算できる．
b) 3 次元角に対する鏡映操作
3 次元角に対しては，2 次元角に対する定義を 3 次元に拡張したものとして考える．すなわち， xnear
i,k∗ を 3 次
元角 k∗ 上の粒子 i との最近傍点として，あとは通常の
near
ポリゴン面に対する場合と同様に計算する．なお，xi,k∗∗は 3 次元角 k 自体の位置に等しい．

- A-04-03 -A-04-03

第26回計算工学講演会(3) 体積補正処理
上述した rwalli,k∗ の定義から，角に対する鏡映操作は，角に関して流体領域を点対称に写像することに等しい．
つまり，角度 θk∗ が 90 度の 2 次元角に対しては，図-5に示すように非整合領域は正しく解消される．しかしながら，θk∗ が 90 度以下になると，角に関して表現している壁領域の体積が非整合領域の体積に対して不足し，90 度以上の場合は反対に体積が過剰に評価される．
角 k∗ の非整合領域を正しく解消するためには，k∗ に関して表現された壁領域の寄与項に補正係数 ϵi,k∗ をかければよい．ϵi,k∗ は，k∗ に関して表現された壁領域のd
体積 Vi,k∗ と，非整合領域の体積 Vi,k
∗ の比として計算できる．
a) 2 次元角に対する体積補正
はじめに，粒子 i が 2 次元角 k∗ 上に存在する場合を∗考える．このとき，|rwalli,k∗ | = 0 を満たすため，ϵi,k は 2 次
元角 k∗ の角度 θk∗ と非整合領域の角度 θkd∗ の比と一致す
る．なお，本稿では角の形状が凹，すなわち 0 ≤ θk∗ ≤ πの場合のみを考える．また，粒子 i の影響半径 h が角∗に接する，すなわち |rwalli,k∗ | = h のとき，Vi,k = 0 かつ
dVi,k
∗ > 0 であるため，ϵi,k∗ は正の無限大に発散する．以上を踏まえると， θd∗(|rwall θkk∗i,k∗ | = 0)(17)ϵi,k∗ = +∞ (|rwall∗ | = h)形領域は三角形に分割できることから，角錐の底面を三角形分割することを考えれば，任意の角推は三角錐に分割でき，その立体角度は式 (21) を用いて計算できる．
式 (18) における平面角度を上述の立体角度に置き換
えれば，3 次元角 k∗ に対する体積補正係数 ϵi,k∗ は次式のように計算できる．

 Ωd
h+η

k∗


 Ωk∗ h+η−|rwall|i,k∗ϵi,k∗ ≈ 


1

(
)
Ω̃i,k∗ < f (Ωk∗ )(22)(else)

ただし，Ωk∗ は 3 次元角 k∗ の立体角度（0 ≤ Ωk∗ ≤ 2π），
Ωdk∗ は k∗ に対する非整合領域の立体角度， f (Ωk∗ ) は式
(18)1 行目の π − θk∗ に対応する Ωk∗ の関数である．ま∗た，Ω̃i,k∗ は rwalli,k∗ と ñk がなす立体角度であり，これは
wall
ri,k∗ を母線，ñk∗ を中心線とした円錐の立体角度として
次式のように計算できる．



∗
rwall


i,k∗ · ñk Ω̃i,k∗ = 2π 1 −
wall|ri,k∗ |(23)Ωdk∗ と f (Ωk∗ ) については，解析的に求めることが困難であるため，本研究では 2 次多項式で近似的に計算する．Ωk∗ = 0 の場合など，特定の条件下で満たすべき値を考慮することで，次式の近似多項式が得られる．

i,k∗と書ける．このとき， rwalli,k∗ の向きが ϵi,k に及ぼす影響を無視すれば，式 (17) から，ϵi,k∗ は次式のように近似
できる．
 θd
h+η

k∗


(θ̃i,k∗ < π − θk∗ ) θk∗ h+η−|rwall|i,k∗ϵi,k∗ ≈ 
(18)


1
(else)Ωdk∗ ≈

258 2Ωk∗ − Ωk∗ + 4π
3π
3

(24)

4 2
11
Ωk∗ − Ωk∗ + 2π
3π
3

(25)

f (Ωk∗ ) ≈ただし，η は式 (10)(11) にもあるゼロ除算を防ぐための安定化パラメータである．θkd∗ は k∗ の角度 θk∗ との幾何的な関係から次式のように計算できる．

(4) 壁面上における境界条件a) 速度のディリクレ条件
壁面上で完全すべり条件を仮定する場合，速度に関する境界条件は以下のディリクレ条件として記述できる．

θkd∗ = 2π − 3θk∗u · n = uwall · n(19)∗θ̃i,k∗ は rwalli,k∗ と ñk がなす角度で，次式のように計算できる．
 wall ri,k∗ · ñk∗ 


(20)
θ̃i,k∗ = 2arccos |rwall|∗i,kただし，ñk∗ は θk∗ を等分する方向の単位ベクトルである．
b) 3 次元角に対する体積補正
上述した 2 次元角に対する体積補正処理を 3 次元に拡張する．まず，3 次元角が有する角度を立体角度として定義する．立体角度の単位はステラジアン (sr) で，0から 4π の値を取る．
三角錐の立体角度 Ω は次式で計算できることが知られている [7]．

Ω=α+β+γ−π© 一般社団法人 日本計算工学会

ここで，uwall は境界面の移動速度である．粒子 i の速度 ui と鏡映粒子 i′ の速度 ui′ の線形内挿で得られる速度分布がこの条件を満たすとすると，鏡映粒子 i′ 上の
速度 ui′ ,k は()
ui′ ,k = Ri,k ui + 2 uwall
·
n
ni,k
i,k
k

(27)

と表せる．
次に，壁面上で非すべり条件を仮定する場合，速度のディリクレ条件は次式のように書ける．

u = uwall(28)

上式を満足するためには，ui′ ,k は(21)

ただし，α, β, γ は錐体を構成する三角形の内角のうち，それぞれ点 O に接するものである．また，任意の多角(26)ui′ ,k = −ui + 2uwallkを満たす必要がある．

- A-04-03 -(29)

A-04-03

第26回計算工学講演会b) 圧力のノイマン条件
圧力に関する境界条件は，壁面上における圧力勾配，
すなわちノイマン条件の形で記述できる．本稿では圧力勾配として静水圧分布を仮定し，次式で表される非一様な圧力のノイマン条件を設定する．

0.10

0.10
0.10

0.20

∂p
= ρg · n∂n

0.10

Unit: m

0.20

(30)

速度の場合と同様に，上式を満たすためには，鏡映粒子 i′ の圧力 pi′ ,k は次式を満たす必要がある．

pi′ ,k = pi − 2ρi g · rwall
i,kPressure evaluation point(31)

Case A(5) 壁面寄与項の計算
以上の諸式を用いれば，ISPH 法の各離散化式における壁寄与項が計算できる．
• 圧力の勾配⟨∇pi ⟩wall =





∑ 
∑


pi 
 p j′ ,kϵi,k ρi Ri,k′ j 
+
∇W



i



2
 ρ2

ρ
j∈Pi′

k∈Wi

j

表–1 解析条件

Time increment, ∆t
Initial particle spacing, l0
Eﬀective radius, hFluid density, ρ0Kinetic viscosity, ν
Relaxing coeﬃcient, αSlip condition on wall surface

(32)

i

⟨
⟩wall
∇2 ui
=


∑ { ( ri′ j · ∇Wi′ j ) (∑  2ϵi,k
)}

ui′ ,k − u j mj
 ρi Ri,k′ j |2 + η2
|r
i
j∈P ′
k∈W
i

(33)
• 速度のラプラシアン (non-slip)
⟨
⟩wall
∇2 ui
=


∑  2ϵi,k ∑ { ( ri′ j · ∇Wi′ j ) (
)}

−
ui′ ,k − u j mj
ρi j∈P ′|ri′ j |2 + η2
k∈W
i

i

(34)
• 圧力のラプラシアン
⟨
⟩wall
∇2 pi
=


∑  2ϵi,k ∑ { ( ri′ j · ∇Wi′ j ) (
)}

mj
pi − p j′ ,k  (35) ρi
′ j |2 + η2
|r
i
j∈P ′
k∈W
i

i

• 速度の発散 (free-slip)⟨∇ · ui ⟩wall =∑  ϵi,k ∑ { (
}
)

m j u j − ui′ ,k · ∇Wi′ j 
 ρi
j∈P ′
k∈W
i

(36)

i

• 速度の発散 (non-slip)⟨∇ · ui ⟩wall =∑  ϵi,k ∑ { (
}
)


m j u j − ui′ ,k · Ri,k ∇Wi′ j 
−
ρ
i
j∈P ′
k∈W
i

© 一般社団法人 日本計算工学会

i

Case B図–6 静水圧問題に用いる容器の形状と寸法• 速度のラプラシアン (free-slip)

i

0.040

(37)

5.0 × 10−4 (s)5.0 × 10−3 (m)3.3l0 (m)3
1.0 × 10 (kg/m2 )1.0 × 10−4 (m2 /s)
5.0 × 10−2
free-slip

ここで，k は角を含めたポリゴン番号を表し，面ポリゴンに対しては ϵi,k = 1 とする．また，Wi は粒子 i の影響半径内に存在するポリゴンの集合である．上式群の導出には，鏡映行列および鏡映粒子の性質として，∇Wi j′ = Ri,k ∇Wi′ j , ρ j′ = ρ j , m j′ = m j ， ri j′ = Ri,k ri′ j ,Ri,k Ri,k = E が成り立つことを利用している．

4. 数値計算例
妥当性検証のため，改良 ERP 壁境界モデルを導入した ISPH 解析コードを用いて 3 次元の静水圧問題を解く．また，比較対象として，角に対する鏡映操作を行わない場合（以降，単に ERP 壁境界モデルと呼ぶ）の計算も実施する．使用する容器の概観を図-6 に示す．容器の形状は角柱で，底面が平らな場合（case A）と，底面に四角錐状のくぼみがある場合（case B）の 2 通りとした．それぞれの図中にある赤点は圧力評価点である．
解析条件は表-1 に示す通りである．

Case A, B における全粒子の水深と圧力値の関係をそれぞれ図-7 と図-8 に示す．図から，case A に対して case
B では，ERP 壁境界モデルを用いた場合の圧力値のばらつきが大きいことが確認できる．これは，底面のポ
リゴンメッシュ形状が複雑になったことにより，壁面寄与の計算誤差が増大したためである．一方，改良 ERP
壁境界モデルを用いた場合，角付近を含めて系全体で滑らかかつ理論値に近い圧力分布が得られている．な
お，case A, B いずれにおいても圧力値が若干過小評価されている粒子が存在しているのは，体積補正処理の近似誤差が原因と思われる．

- A-04-03 -A-04-03

第26回計算工学講演会

1200

ERP
Improved ERP
Theory3200
Pressure (Pa)

1000
Pressure (Pa)

4000

ERP
Improved ERP
Theory

800
600
400
200

2400
1600
800

0
0

0.02

0.04 0.06
Depth (m)

0.08

0

0.1

図–7 水深と圧力の関係（case A）

1600
Pressure (Pa)400
00.02 0.04 0.06 0.08
Depth (m)

0.1

0.7
0.8
Time (s)0.91

操作を定義し，近似的な代数計算による壁領域の体積補正処理を導入した改良 ERP 壁境界モデルを開発した．
その際，角を 2 次元角と 3 次元角に分類し，2 次元角に対して定義した鏡映操作および体積補正を立体角度の概念を用いて 3 次元角に拡張した．
続いて，提案手法を ISPH 法に導入し，数値解析例として 3 次元の静水圧問題を解いた．計算結果から，従来の ERP 壁境界モデルを用いた場合は角付近で不自然な圧力撹乱が生じたのに対し，改良 ERP 壁境界モデルを用いた場合は圧力撹乱が解消され，系全体で滑らかかつ理論値に近い圧力分布が得られた．

800

0

0.6

図–10 容器底面の評価点上における圧力時刻歴（case B）

ERP
Improved ERP
Theory

1200

0.5

0.12

図–8 水深と圧力の関係（case B）

参考文献
2000
ERP
Improved ERP
TheoryPressure (Pa)

1600
1200
800
400
0

0.5

0.6

0.7
0.8
Time (s)0.91

図–9 容器底面の評価点上における圧力時刻歴（case A）

次に，case A, B における評価点上の圧力時刻歴を図-9と図-10 に示す．図-9, 10 から，case A, B いずれにおいても，ERP 壁境界モデルを用いた場合に時間的な圧力
撹乱が生じていることがわかる．一方，改良 ERP 壁境界モデルではいずれのケースにおいても安定的に圧力を計算できている．なお，計算結果が理論値を全体的に下回っているのは，密度一定条件の緩和係数 α により，微小に圧縮性を許容していることが原因である．

5. 結論
本研究では，従来の ERP 壁境界モデルにおいて生じる角付近のエラーを解消するために，角に対する鏡映© 一般社団法人 日本計算工学会[1] Monaghan, J. J., Smoothed particle hydrodynamics,
Journal of Computational Physics, Vol.110, No.1,
pp.543-574, 1992.
[2] Koshizuka, S. and Oka, Y., Moving-particle semiimplicit method for fragmentation of incompressible ﬂuid, Nuclear Science and Engineering, Vol.123,
pp.421-434, 1996.
[3] Attaway, S., Heinstein, M. and Swegle, J., Coupling ofsmooth particle hydrodynamics with the ﬁnite element
method, Nuclear Engineering and Design, Vol.150,
No.2, pp.199-205, 1994.
[4] Mitsume, N., Yoshimura, S., Murotani, K. and Yamada, T, Improved MPS-FE ﬂuidstructure interaction
coupled method with MPS polygon wall boundarymodel, Computer Modeling in Engineering and Sciences, Vol.4, pp.229-247, 2014.
[5] Mitsume, N., Yoshimura, S., Murotani, K. and Yamada, T., Explicitly represented polygon wall boundary model for the explicit MPS method, Computational Particle Mechanics, Vol.2, No.1, pp.73-89,
2015.
[6] Asai, M., Aly, A. M., Sonoda, Y. and Sakai, Y., A
stabilized incompressible SPH method by relaxing thedensity invariance condition, Journal of Applied Mathematics, Vol.2012, 139583, 2012.
[7] Frankl, P., 前原潤, 幾何学の散歩道 - 離散・組み合わせ幾何入門, 共立出版, 1991.

- A-04-03 -