Vol.29 (2024年6月)

計算工学講演会論文集

日本計算工学会

粒子法とASI-Gauss法を用いた
非解像型カップリング手法の開発
Development of an Unresolved Fluid-Structure Coupling Scheme using Particle Method
and ASI-Gauss Method
大日向智志1)，大村浩之2)，浅井光輝3)，磯部大吾郎4)
Satoshi Ohinata, Hiroyuki Omura, Mitsuteru Asai and Daigoro Isobe
1) 筑波大学大学院（〒305-8573 茨城県つくば市天王台1-1-1，E-mail: s2420854@u.tsukuba.ac.jp）
2) 博（工）防災科学技術研究所 特別研究員（〒305-0006 茨城県つくば市天王台3-1, E-mail: homura@bosai.go.jp）
3)博（工）九州大学 工学研究院 准教授（〒819-0395 福岡県福岡市西区元岡744, E-mail:asai@doc.kyushu-u.ac.jp ）4)博（工）筑波大学 システム情報系 教授（〒305-8573 茨城県つくば市天王台1-1-1, E-mail:isobe@kz.tsukuba.ac.jp）
In this paper, an unresolved fluid-structure coupling scheme was developed. The ISPH method, which isa particle method for incompressible flows, was used as a numerical method for fluid analysis. The ASIGauss method, which can stably calculate nonlinear phenomenon including fracture was used as anumerical method for structural analysis. The numerical results had shown that the proposed method cancalculate fluid forces reasonably well even when particle resolution is coarse with respect to member
dimension. Furthermore, the proposed method can simulate process of fracture and advection of beam
member subjected to fluid forces.
Key Words : ISPH method, ASI-Gauss method, FSI Problems, Unresolved coupling scheme

1． はじめに
日本は世界有数の地震大国であり，2011年3月11日に発生した東北地方太平洋沖地震など，これまでに多くの津行う．流体解析には粒子法のうち，非圧縮性流れに適用できるISPH法[5]を用いる．構造部材の解析には，破断の表
現を精度よく行えるASI-Gauss法[6]を用いる．

波被害を受けている．津波対策を講じる際，津波のような自由表面を伴う非圧縮性流れを解く数値解析手法として

2． ISPH法

粒子法が有用である．さらに，津波被害をより正確に予測(1) 基礎方程式

するためには，建物の流失など，流体によって構造が破壊，

非圧縮性流れの質量保存則およびラグランジュ型の運大移動する過程まで考慮する必要がある．そのため，津波

動量保存則は以下の式のように書ける．

と建物が相互に影響しあう流体-構造連成 (Fluid-StructureInteraction，FSI) 問題として扱う必要がある．粒子法を用いたFSI解析手法はこれまでにいくつか提案されている([1][2]など) が，いずれも構造物の形状を詳細に表現する，
つまり粒子で構造物を空間的に解像することを前提としている (以下，このようなFSI解析手法を解像型と呼ぶ) ．
一方，都市レベルで浸水域を予測するような大規模な
津波遡上解析では，現実的に使用できる計算資源を踏まえると粒子の直径を数[m]程度の大きさにすることが一般的である[3][4]．この粒子径オーダーに対して，樹木や

∇･𝒗
1
∇𝑝
𝜌

𝐷𝒗
𝐷𝑡

0

1

𝜈∇ 𝒗

𝒈

2

ここで，𝒗は速度，𝑡は時間，𝜌は密度，𝑝は圧力，𝜈は動粘性係数，𝒈は重力加速度である．
(2) 時間積分スキーム
ISPH 法では，速度を陽的に，圧力を陰的に計算する射影法によって各ステップでの速度，圧力を計算する．
速度の更新は以下の式により行う．

家屋の構造部材などの断面寸法は高々数十[cm]程度であるため，部材の形状を十分細かく解像することができな
い．したがって，既存の連成解析手法はこのような広域の問題には適用できない．
そこで本研究では，部材の形状を陽に表現せずに，部材

𝒗∗
𝒗

𝒗
𝒗∗

𝜈𝛻 𝒗

𝒈 Δ𝑡

3

1
𝛻𝑝
𝜌

Δ𝑡

4

に作用する流体力を抗力モデルに基づき流速から計算すここで，添え字∗は n ステップと n+1 ステップの間に定義る非解像型のFSI解析手法を開発し，その妥当性の検証をされる中間状態を表す．また，Δ𝑡は時間増分である．

圧力の更新は次式の圧力ポアソン方程式を解くことで行う．

とる．この特性を利用し，次式に従って自由表面判定を行
う．

𝜌
𝛻･𝒗∗
𝛥𝑡

𝛻 𝑝

𝛼 𝜌 𝜌
𝛥𝑡

𝜆
𝜆

5

ここで，𝜌は粒子の数値的な密度，𝛼 は緩和係数である．

𝜆

→𝑖∈𝔽
→𝑖∉𝔽

𝜆

11

ここで，𝔽は自由表面に属する粒子の集合を表す．また，(3) SPH 法に基づく空間離散化

本研究では𝜆領域Ω内で連続なスカラー場𝜑は関数のカーネル近似

ずれの場合も満たさない粒子に対しては，次の段階の幾により次式のように近似できる．

0.9とした．式 (11)のい何的判定に移る．
幾何的判定では，式 (12) ，(13) で定義した外向き法線

⬚

𝜑 𝒙

0.2，𝜆

𝜑 𝝃 𝑊 𝝃𝒙 𝑑𝝃6ベクトルの方向の近傍に定義する探索領域に別の粒子が存在しない場合，自由表面に属する粒子であると判定す

SPH 法はこのカーネル近似に基づいた離散化スキームでる．

あり，式 (6) を粒子離散化近似した形で書くと，
𝒏
〈𝜑〉

𝜑 𝒙

𝑉𝜑 𝑊

〉

12

〉

7

∈ℙ

〈𝛻𝜆

〉

𝒎

𝑉 𝜆

𝜆 𝛻𝑊

13

∈ℙ

となる．ここで，ℙ は粒子𝑖の影響半径内に存在する粒子の集合，𝑉 は粒子𝑗が代表する体積である．

(5) 壁面の表現

勾配，ラプラシアンなどの微分表現は，式 (7) と𝜑のテイラー展開から導出され，それぞれ以下の式で表される．
勾配は，

本研究では水路の床，壁面の表現には大村ら (2023) [8]の IERP 壁境界モデルを用いる．IERP モデルは，陽解法の粒子法への適用が前提の ERP モデルをベースに，角部分の表現を精緻化し，半陰解法の粒子法に適用できるよ

〈𝛻𝜑〉

𝛻𝜑 𝒙

𝑉 𝜑

𝜑 𝛻𝑊

8

∈ℙ

(1) 離散化スキーム
空間の離散化はすべての変位が線形関数で定義される

〈𝛻 𝜑〉
2

うに発展した壁境界モデルである．

3． ASI-Gauss法と書ける．同様にラプラシアンは，
𝛻 𝜑 𝒙

〈𝛻𝜆
〈𝛻𝜆

𝑉 𝜑

𝜑

∈ℙ

𝒙

𝒙

𝒙

𝒙

𝛻𝑊

9

線形チモシェンコはり要素を用いる．また，時間積分法にはNewmarkのβ法を用いる．これらを用いると運動方程式，
速度，加速度ベクトルは以下の式のように与えられる．

と書ける．

𝑴𝒖

(4) 自由表面判定

𝑪𝒖Δ

自由表面での圧力の境界条件を付与するため，自由表
𝒖面に位置する粒子を判定する必要がある．本研究では自

1

Δ

Δ

1
𝒖
2𝛽

由表面判定に Marrone ら (2010) [7]の方法を用いる．この方法では，モーメント行列𝒎 と呼ばれる周囲の粒子分布を表す行列およびその最小固有値𝜆を用いた判定と，

𝒖

Δ

1

𝛿
𝒖 Δ𝑡
2𝛽

𝑲Δ𝒖
1
𝛽Δ𝑡
1

𝑭

𝒖

Δ

1
𝛽Δ𝑡

𝛿
𝒖
𝛽

𝑹

14

Δ𝒖

15

𝛿
𝛽Δ𝑡

Δ𝒖

16

幾何的な判定の 2 段階で判定を行う．
モーメント行列𝒎 は次式のように表される．

ここで，𝑴は全体質量マトリックス，𝑪は全体減衰マトリックス，𝑲は全体剛性マトリックス，𝑭は節点外力ベクト𝒎

𝑉 𝒙𝒙 ⊗ 𝛻𝑊

10∈ℙ

周囲に粒子が規則正しく充満している場合，𝒎 は単位行列となり，𝒎 の最小固有値は𝜆

ル，𝑹は節点内力ベクトルである．また，𝛽，𝛿はNewmarkのβ法における重みパラメータであり，本研究ではLynnら

1となる．一方，自由表面付近など，粒子が疎な位置では𝜆は 0 に近い値を(2007) [6]の研究を参考に𝛽

4/9，𝛿

5/6とする．

ASI-Gauss法では，1つの梁部材を2つの線形チモシェン
コはり要素で表現する．弾性域では応力評価点が3次はり要素のものと一致するように数値積分点を配置し，塑性

ヒンジが生じた後は塑性ヒンジの位置と応力評価点が一の作用点はそれぞれの小部材の重心位置とした．

致するように数値積分点をシフトすることで，曲げ変形d)

小部材ごとの流体力計算

式 (18) の抗力の式を小部材𝑖に働く抗力の式に書き換を精度よく解ける．
(2) 梁部材接触モデル

えると次の式となる．

梁部材同士および梁部材と壁面の接触の表現は文献[8]の粒子ベース接触モデルを利用する．粒子ベース接触

1
𝐶 𝜌𝒗 𝑐𝑜𝑠𝜃|𝒗 𝑐𝑜𝑠𝜃|𝐷𝑙
2𝒇

モデルでは部材形状を粒子によって表すため，部材の形状を考慮した接触が表現できる．

𝒇 は小部材𝑖に働く抗力，𝒗 は小部材𝑖における流速，𝑙 は(3) 梁部材破断モデル
破断の判定は文献 [8]の研究を参考に，以下の式で表されるひずみ量による判定式を用いる．
𝜒
𝜒

𝜒
𝜒

𝜀
𝜀

𝜀
𝜀

小部材𝑖の水没部長さを表す．
小部材𝑖における𝒗 ，𝜆を設置し，計算点𝑖で内挿近似を行う．本研究では，計算点𝑖

𝜀
𝜀

1

17

ん断ひずみである．また，上添え字の𝑐はひずみの破断臨

小部材𝑖の水没部長さ𝑙 は，計算点𝑖(1) 抗力モデルで計算された𝜆の値を利用し，以下の式に従って計算する．

𝑙

⎧
⎪ 𝜆
𝑑
⎨ 𝜆
⎪
⎩

にを破断面とし，新たな節点を生成して要素を分離させ

4． ISPH法とASI-Gauss法の非解像型カップリン
グ手法

での内挿計算には自由表面や壁近傍の粒子の空隙域を考慮し，修正カーネル関数を利用する．

界値を表す．破断が判定された要素は，塑性ヒンジの位置
る．

の設置位置

を小部材の重心位置から部材幅分上流の位置とした．計
算点𝑖率，𝜀 は軸ひずみ，𝜀 ，𝜀 はそれぞれ𝑥軸，𝑦軸方向のせ

0

𝜆

𝜆
𝜆
𝑑

𝜆

𝜆
𝜆

𝜆

𝜆

20

𝜆

ここで𝑑は小部材の長さである．
抗力係数に関して，本研究では対象とする梁部材を正梁部材に作用する流体抗力

方形断面の角柱と考え，Sevhan ら (2018) [10]の研究を参山邊ら (2005) [9]を参考に，流れ方向と梁部材のなす角𝜃を用いると，流れの中に存在する梁部材に働く抗力𝒇 を
以下の式 (18) のように計算する．
𝒇を計算する際は，小部材𝑖の重心位置から上流側に計算点𝑖

ここで，𝜒 ，𝜒 はそれぞれ要素座標系𝑥軸，𝑦軸方向の曲a)19

考に角柱の迎え角に応じた抗力係数の近似曲線を用いた（図 1）．
また，小部材に働く流体力として，抗力の他に浮力を考1
𝐶 𝜌𝒗cos𝜃|𝒗cos𝜃|𝐷𝑙2

える．浮力は以下の式により計算する．
18
𝒇

𝜌𝒈𝐴𝑙

21

ここで，𝐶 は抗力係数，𝒗は流速，𝐷は梁部材の流れ方向
への投影幅，𝑙は梁部材の水没部長さである．
b)

ここで𝐴は部材の断面積，𝑙 は小部材𝑖の水没部長さ，𝒈は
重力加速度である．

水没部長さ式 (18) を用いて抗力を計算するためには，梁部材の水構造解析の手法に FEM を利用しているため，以上によ没部長さ𝑙を求める必要がある．そこで，式 (10) のモーメ

り求めた小部材ごとの流体力𝒇 ，𝒇 をそれぞれの要素のント行列𝒎の最小固有値𝜆を利用し，部材の水没部長さを計算する．平らな自由表面上での𝜆を考えると，影響半径の半分は水粒子で満たされ，もう半分は水粒子が存在しない空隙域となる．その場合，𝜆は理論上𝜆0.5となる．したがって本研究では，梁部材のうち𝜆𝜆(𝜆0.5 ) となる部分を水没しているとみなし，部材の水没部長さ𝑙を計算する．
c)

部材の分割

はり部材の軸方向に流速が分布する場合，抗力の値も流速分布に応じた分布を持つ．そこで，はり部材を流体の
空間解像度に応じて複数の小部材に分割し，それぞれの小部材ごとに流体力計算を行う．本研究では，分割間隔を流体粒子の粒子径𝑑とした．また，小部材が受ける流体力

図1

迎え角-抗力係数関係の近似曲線

節点力に変換する必要がある．小部材𝑖の抗力𝒇 の節点力確認できた幅0.05[m]のものと比較することで，角柱に対変換は以下の式 (22) ，(23) に従って行う．また，浮力𝒇する粒子の解像度ごとの妥当性の検証を行う．なお，変形も同様に変換する．

の影響を無視するため，角柱の断面定数は前項で使用した0.05[m]の角柱のものを用いる．

𝑭

𝑭

1
1
2

𝜉 𝒇

1
1
2

𝜉 𝒇

解析の結果を図5に示す．図5より，正規化した抗力が角22

柱幅0.05[m]の場合と近い値をとっていることが確認でき23

る．流体力が発生するタイミングに遅れが生じている原因は，抗力計算に用いる物理量を計算する計算点𝑖の位置の違いによるものである．すなわち，計算点𝑖はここで添え字 1，2 は局所節点番号を表す．また，𝜉 は小小部材の重心位置から上流側に角柱幅分シフトさせた位部材の重心位置を要素座標系で表したもので， 1

置に設けているため，角柱幅が大きいケースほど早いタ𝜉

1の値をとる．
e)

イミングで流体力が発生しているためである．また，0.2[s]流体への反作用力付近の抗力のピーク値が部材幅0.05，0.025[m]のケースで流体と構造物間の力のつり合いを満たすため，𝒇 の逆は部材が細くなるにつれ抗力が小さくなっている一方，向きの力 𝒇 を小部材𝑖の重心を中心とする影響範囲内

部材幅0.01~0.0005[m]のケースでは部材幅が細くなるにとして加える．な

つれて抗力が大きくなっている．前者の理由は計算点の流体粒子𝑘へ分配して，反作用力𝒇
お，反作用力𝒇𝒇は以下の式により定義する．
𝒇

𝑉𝑊∑ ∈ℙ 𝑉 𝑊𝑖の位置は小部材の重心位置から部材幅分シフトさせた位置に設けいているため，水没部長さ𝑙の値が部材幅𝑘∈ℙ24が大きいほど大きくなったからである．一方後者の理由
は，部材から流体への反作用力により部材付近の流速が低下していることが原因である．部材幅が粒子径に対し5． 提案手法の妥当性検証て大きいほど大きな反作用力が流体粒子に作用し流速の(1) 解像型連成解析との比較低下も大きくなる．部材幅0.05，0.025[m]のケースでは部

提案手法の妥当性を検証するため角柱部材付きダムブ

材幅が十分大きいため，部材付近の流速低下の影響を受レイク問題を，提案手法である非解像型連成手法とPMS

けない位置に計算点𝑖モデル[11]による解像型連成手法それぞれで解き，角柱に0.01~0.0005[m]のケースでは部材幅が小さく計算点𝑖

働く流体力を比較する．

を部材付近の流速低下の影響を受けない程度上流側にシを設置できていたが，部材幅対象とする系の概観を図2，3に示す．流体の粒子径はフトできていなかった．実際，小部材ごとの抗力計算に用10 [m]で固定とした．この粒子径は一辺0.05[m]の角いた流 速の 平 均値の 時刻 歴 (図6) を見ると，部 材幅

柱の断面形状を十分詳細に表現できるため，解像型連成0.01~0.0005[m]では部材幅が細くなるごとに平均流速が手法により精度よく角柱に働く流体力が計算できる．角大きくなることが確認できる．この結果から，多少の差が柱の材質は鋼材を想定し材料定数を設定した．流体解析

存在するが，粒子径が角柱の形状に比べて粗い場合にお条件，構造解析条件はそれぞれ表1，2の通りである．水路いても妥当に流体力計算が行えることが確認できた．

5の壁面と流体の間の摩擦は考慮しない自由滑り条件とした．また，ダムブレイク問題では流れの方向がX軸正の方向であるため，非解像型連成手法ではX軸負の方向を上流側の方向として物理量計算点を設置した．
解像型連成手法と提案した非解像型連成手法の結果を比較したグラフを図4にまとめた．図4をみると，提案手法
での結果は解像型連成手法の結果によく一致しており，

図2

ダムブレイク問題における系の概観(正面図)

図3

ダムブレイク問題における系の概観(側面図)

提案手法により妥当に流体力が得られることが確認できる．
(2) 空間解像度が異なるケースでの比較
部材の断面寸法に対する流体の空間解像度の比が異なるケースの比較を行う．検証対象の系は5. 1節と同様のものとし，粒子径は固定のまま角柱幅を変える．角柱の一辺の長さは0.0005，0.001，0.0025，0.005，0.01，0.025，0.05[m]の7種類とする．得られた流体力を角柱幅𝐷で除すことで正規化する．正規化したそれぞれの値を前項で妥当性が表1

ダムブレイク問題における流体解析条件
5.0

時間増分

10
1000

総ステップ数
5.0

粒子径

s
m

10

128000

総流体粒子数
影響半径

粒子径

緩和係数

5.0

2.4

10 
1.0

滑り率
密度

1.0

10

kg/m 

動粘性係数

1.0

10

m /s 

9.8 m/s 

重力加速度

表2

ダムブレイク問題における構造解析条件

時間増分

5.0

10

s
N/m 

7.93

10

降伏応力

5.98

10

N/m 

密度

4.91

10

kg/m 

0.4

ポアソン比

重力加速度

抗力計算に用いた X 軸方向流速の平均値の角柱幅ごとの比較

1000

総ステップ数ヤング率

図6

9.8 m/s (a) 時刻 0.175[s]で角柱が破断する様子

図4(b) 時刻0.300[s]で角柱が移流する様子

PMS モデルと提案手法の X 軸方向流体力の比較

図 7 角柱流失シミュレーションの様子

6． 角柱の流失シミュレーション
提案手法により角柱の流失シミュレーションを行う．
対象とする系は5.1節で用いたものと同様とする．また，角柱の破断が生じるように構造部材の降伏応力の値を仮想的に小さくする．さらに，角柱が破断後に流れに沿って移流する過程を表現するために，部材軸方向の抗力も計算する．簡単のため，軸方向の抗力計算に用いる抗力係数は軸直交方向の抗力係数と同じものを使用した．シミュレーションの可視化結果を図7に示す．図7を見ると，角柱が流体により破断し，流されていく様子が再現されている．
図5

正規化した X 軸方向流体力の比較

7． 結論domaindecomposition

explicit

MPS

methodand

本研究では，粒子径が部材の寸法に対し大きい場合にapplication to large-scale tsunami analysis with floating

おいても発生する流体力を妥当に計算できる非解像型のobjects, Journal of Advanced Simulation in Science and

流体-構造連成解析手法の開発を行った.妥当性および有効性を検証するために，開発した非解像型の連成手法と，

Engineering, Vol. 1, No.1, pp. 16-35, 2014.
[5] M．Ellero, M. Serrano, and P. Español, Incompressible

PMSを利用した解像型の連成手法をダムブレイク問題にsmoothed

より比較したところ，以下の結果が得られた．

Computational Physics, Vol. 226, No. 2, pp. 1731-1752,●●

解像型連成では計算が不可能なほど部材断面寸法

particlehydrodynamics,

Journalof

2007.

に対して粒子径が粗い場合においても，提案した[6] K. M. Lynn and D. Isobe, Finite Element Code for Impact

非解像型連成では妥当に流体力を計算できること

Collapse Problems of Framed Structures, Internationalが確認された．

Journal for Numerical Methods in Engineering, Vol. 69,

提案した非解像型連成解析手法によって，構造部

No. 12, pp. 2538-2563, 2007.

材が流体力を受け，破断，流出していく挙動を再現できることが確認された．

[7] S. Marrone, A. Colagrossi, D. Le Touzé, and G. Graziani,
Fast free-surface detection and level-set function definition
in SPH solvers, Journal of Computational Physics, Vol.

参考文献

229, No. 10, pp. 3652-3663, 2010.

[1] C. Antoci, M. Gallati and S. Sibilla, Numerical simulation[8] 大村 浩之，三目 直登，浅井 光輝，磯部大吾郎，粒of fluid-structure interaction by SPH, Computers &子法とASI-Gauss法のカップリングに基づく波浪-骨Structures, Vol. 85, pp. 879-890, 2007.

組構造連成解析，計算工学講演会論文集，Vol. 28,[2] G. Fourey, G. Oger, D. Le Touzé and B. Alessandrini,

2023.

Violent fluid-structure interaction simulations using a[9] 山邉 茂之，熊本 博光，西原 修，釣り仕掛けアニメ

coupled SPH/FEM method, IOP Conference Series:

ーションのための抗力係数の同定と過渡的糸挙動のMaterials Science and Engineering, Vol. 10, p. 012041,

再現，精密工学会論文集，71巻，11号，pp. 1448-1453,

2010.

2005.

[3] 江口 史門，浅井 光輝，大谷 英之，一色 正晴，建物[10] M. Sevhan, M. Sarioglu, Y. E. Akansu, Effect of Attack群を含む地表面詳細モデルを用いた粒子法による三Angle on Flow around A Square Prism with A Splitter

次元津波遡上解析，土木学会論文集A1（構造・地震

Plate, Uludağ University Journal of the Faculty of工学），Vol. 72, No. 4（地震工学論文集第35巻）, I_367-

Engineering, Vol. 23, No. 2, 2018.

I_377, 2016.

[11] S. Koshizuka, A. Nobe, and Y. Oka, Numerical analysis of[4] K. Murotani, S. Koshizuka, T. Tamai, K. Shibata, N.

breaking waves using the moving particle semi-implicit

Mitsume, S. Yoshimura, S. Tanaka, K. Hasegawa, E.

method, International Journal for Numerical Methods in

Nagai, and T. Fujisawa, Development of hierarchical

Fluids, Vol. 26, No. 7, pp. 751-769, 1998.