第 24 回数値流体力学シンポジウム
F10-4

3 次元 SPH 法による非圧縮性流体と弾性体の連成シミュレーション
Three-dimensional SPH Simulation of Incompressible Fluid Interacting with an Elastic Structure⃝

佐竹 正俊, 東北大・院, 仙台市青葉区荒巻字青葉 6-6-01, E-mail：satake@ﬂuid.mech.tohoku.ac.jp
茂田 正哉，東北大・工, 同上, E-mail：shigeta@ﬂuid.mech.tohoku.ac.jp
伊澤 精一郎，東北大・工, 同上, E-mail：izawa@ﬂuid.mech.tohoku.ac.jp
福西 祐，東北大・工, 同上, E-mail：fushi@ﬂuid.mech.tohoku.ac.jp

Masatoshi SATAKE, Tohoku University, Aramaki-Aoba 6-6-01, Sendai, JapanMasaya SHIGETA, Tohoku UniversitySeiichiro IZAWA, Tohoku University
Yu FUKUNISHI, Tohoku University
Fluid-structure interaction model for the three-dimensional SPH method is proposed and applied to a liquidcolumn-break problem with an elastic gate. First, the oscillation of a cantilever is simulated using an elastic
analysis model, which is extended from the two-dimensional conventional model to a three-dimensional one. The
frequencies close to the theoretical values are obtained in the simulation. Next, the fluid-structure interaction
model is made by coupling the elastic analysis model to a model analyzing an incompressible flow. The force
exerted against the structure by the fluid is obtained using the traveling distances of the water particles during
the density adjusting procedure. The computational result using the new model agrees well with the results ofexisting experiments and simulations.

1.

緒言

計算方法
SPH 法においては密度は以下のように周囲粒子の質量の重ね合わせとする．
∑ρ(ri ) =mb Wab(1)

2.

タイヤの表面にある溝や切り込み（トレッドパターン）
は，運動性能や静粛性，燃費性能といったタイヤの基本性能を決める重要な要素である．例えば，濡れた路面における排水性能を向上させるには太い溝が適しているが，ロードノイズを低く抑えるためには溝は細い方が望ましい．このように相反する要求を満たす製品を開発するた
めには，現象を正確に捉え，運動を支配するメカニズムを理解することが必要不可欠である．しかし，タイヤが水の溜まった路面を走行するような場合，タイヤの通過にともなって水面は大きく変化し，またトレッドに沿って巻き上げられた水が排水されるなど，現象は極めて複雑である．さらにタイヤは自重により変形するので，剛体として扱うことはできない．粒子法の 1 つである SPH（Smoothed Particle Hydrodynamics）法は，このような複雑な流れ場であっても，適切な物理モデルを組み込むことで適用できる可能性がある．

SPH 法はもともと天体物理学の分野で発達した計算手法で，密度を媒介変数とすることから，原理的には弱圧縮性流体の数値解法である．アルゴリズムが簡便で固体も流体同様の扱いができることから，構造体の変形や破壊といった問題にも適用が進み，流体工学の分野のみならず土木や材料関係の分野でも幅広く使われている．例
えば，関根ら 1) は積層板の解析モデルを構築し，非接着で重ね合わせた単純な積層構成のアルミニウム積層板に対して，超高速で飛翔するアルミニウムの飛行体が衝突した際の変形破壊挙動について計算を行っている．一方
我々のグループでは，これまでに非圧縮流体に拡張したSPH 法 2) を二重拡散対流 3) や液柱の崩壊 4) ，スロッシ
ング 4) ，液滴の衝突 5) や相変化のある流れ 6) などへの適用を進めている．本研究では，これらの流体解析モデルに既存の弾性体解析モデルを連成させることで，SPH法を用いて濡れた路面を走行するタイヤのような複雑な流体挙動を示す現象を解析することを目指している．本稿では，まず Gray ら 7) の弾性体解析モデルを 3 次元に拡張して片持ち梁の振動解析を行い計算結果の妥当性について検証を行った後，弾性板のゲートを通過する液柱の崩壊現象を取り上げ，非圧縮流体解析モデルと連成させた計算を行ない，現象の再現を試みた．

b

また，運動方程式は以下のようになる．
()∑pbpa ∂Wabdvaimb=−+ 2dtρ2bρa∂xjab(2)

ここで a, b は粒子番号，ρ は密度，m は質量であり，Wはカーネル関数を表す．本計算では，カーネル関数として一般的な M4 スプライン関数を用いる．

2.1 弾性体解析
式 (1) をそのまま用いると，外部との境界面近傍で密度が急激に低下し非現実的な圧力勾配が生じる．そこで，密度の異なる流体との連成を視野に入れて，次式のように連続の式を微分形式に書き換える．
∑ mbdρai ∂Wab= ρav̂ba(3)dtρb∂xiabv̂ba = v̂a − v̂b は粒子 a 粒子 b 間の平均化された相対速度の差であり，平均化された相対速度は次式で計算される．
∑mcv̂ai = vai + ϵ̃(vci − vai )Wac(4)1(ρ+ρ)
acc 2ϵ̃ (0 ≤ ϵ̃ ≤ 1) は粒子 a の速度をその周囲にある粒子 c の速度を用いてどの程度平均化するかを決めるパラメータであり，Gray ら 7) と同様に ϵ̃ = 0.5 とした．なお，式(3) の微分形式表示では，質量が厳密には保存されなくなることに注意が必要である．
弾性体粒子に対する運動方程式は, Gray ら 7) の方法に従い外力項を消去した式 (2) の運動方程式に人工粘性項と人工応力項を付加したものを用いる．
()∑σbijdvai
σaij∂Wabij nij=mb(5)+ 2 + Πab δ + Rab f2dt
ρaρb∂xjab1c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
F10-4

ここで σaij ，σbij は粒子 a，b の持つ応力テンソルであり，ij nΠab δ ij ，Rabf はそれぞれ人工粘性項，人工応力項である．人工粘性は粒子同士の異常接近を抑えるために，人工応力は Tensile Instability を排除するために導入されたものである．

2.2 非圧縮性流体解析
式 (1), 式 (2) を基礎方程式とする SPH 法は，密度を媒介変数とする圧縮性流体の数値解法である．非圧縮性流体は密度変化を許容しない”硬い”流体であるから，このまま式を適用することはできない．本研究では，以下の手順に従い，場の密度変化がある一定レベル以下に収まるように粒子を運動させることで近似的に非圧縮性の実現を図る 2) ．
まず，各粒子を慣性に従って Euler 的に移動させる．
r∗ia (t + ∆t) = rai (t) + vai (t)∆t

3.

(6)

(7)

簡便のため基準圧力 p0 をゼロとして粒子の圧力 pa を求め， 粒子に加わる圧力勾配を計算し，(2) 式の運動方程式に従って粒子位置の補正を繰り返す．粒子密度の変動が場の平均密度の 1% 以内に収まった段階で非圧縮性が達成されたものとみなし，このときの粒子位置を次時刻の粒子位置 rai (t + ∆t) とする．移動距離より速度を以下のように求めv ∗ia (t + ∆t) =rai (t + ∆t) − rai (t)∆tva2 =))Vf ( (M cos(kx1 ) − cosh(kx1 ) .
Q( ())− N sin(kx1 ) − sinh(kx1 )

(12)

ここで，

M = sin(kL) + sinh(kL)(8)

N = cos(kL) + cosh(kL)

さらに粘性項を加えて，次時刻の速度は以下のように定める．

vai (t + ∆t)

問題設定

3.1 片持ち梁の振動解析
弾性体の運動解析の対象として，片持ち梁の振動を取り上げた．梁は振動が起こりやすいように，ヤング率を低めに設定して計算を行った．Fig. 1 に初期粒子配置を示す．緑色で着色した粒子は静止した剛体壁を，灰色で着色した粒子は弾性体の梁を構成する粒子を示している．計算には，壁面粒子 336 個，弾性体粒子 960 個の合計 1,296個の粒子を用いた．梁の寸法は，長さ 40mm (=L)，高さ4mm (=H)，幅 6mm で，粒子の直径は 1.0mm である．
また，弾性体のヤング率 E は 1.859MPa，ポアソン比 νは 0.49 とした．Gray らと同様，初期条件として +x2 方向に次式で示される速度 va2 を与えた．

ここで添え字の ∗ は，仮の位置であることを示す．粒子は単に現在の速度で移動しただけなので，粒子の分布には偏りが生まれ局所的な密度勾配が発生する．等エントロピー変化を仮定すると，密度の変化 ∂ρ は音速 c を介して圧力の変化 ∂p と関連付けることができる．
)(∂ρa= c2
∂pa s∴ pa = p0 + c2 (ρa − ρ0 )そして，以上の手順で求めた流体から弾性体に作用する単位質量あたりの力（加速度）αai f →s を，弾性体の運動方程式の外力項に加えることで，両者の運動を連立させた．なお係数 γ の値は γ = 10.0 とした．
また，流体・弾性体の連成解析では，時間進行についても注意が必要である．弾性体解析では圧縮性，流体解析では非圧縮性を前提に計算を行うため，両者の時間刻みは大きく異なる．例えば，流体解析では，密度補正プロセスを経ることで時間刻みを大きくとることができるが，弾性体解析では，その 200 分の 1 にまで時間刻みを落とさなければ安定した計算を行うことはできない．そこで本計算では，流体解析の時間刻みを弾性体解析の 200倍とし，200 回に 1 回両者を同時に解く方法を採用した．
すなわち，流体と弾性体の間に作用する力は，200 回ごとのみに更新される．

Q = 2 (cos(kL) sinh(kL) − sin(kL) cosh(kL))= v ∗ia (t + ∆t)i∑(µa + µb )vab∂Wab+mb∆t (9)iρρr∂xjaababb(13)であり，梁の自由端（x1 = L）における速度 Vf は 0.01とした．また，k は，式cos(kL) cosh(kL) = −1

2.3 流体・弾性体の連成解析
流体と弾性体の運動を連成させる際にポイントとなるのは，力の取り扱いである．本計算では，非圧縮性流体解析において粒子位置を補正する際に加えられた局所的な密度勾配に起因する力を積分し，流体と弾性体に作用する力とした．具体的な計算手順は以下の通りである．

(14)を満たす値である．
Rigid wall
Elastic beam

(1) 密度補正プロセスを通して水粒子 a の移動した距離から単位質量あたりの力（加速度）άai を算出
irai (t + ∆t) − ra∗ (t)iάa =2(∆t)(10)(2) 弾性体表面粒子 b の法線ベクトル nib に対する άbi の正射影ベクトル αbi を計算し，カーネル関数による足し合わせを行う.
∑ mb
αai f →s ＝γάi Wab (−nia )
(11)
ρb b

Fig. 1: Initial particle distribution of the cantilever
problem.

b2c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
F10-43.2 弾性板のゲートを通過する液柱の崩壊現象
流体と弾性体の連成問題としては，弾性板のゲートを通過する液柱の崩壊現象を取り上げた．Figure 2 に粒子の初期配置を示す．弾性体解析の場合と同様に緑色の粒子は壁粒子を灰色の粒子は弾性体粒子を表しており，濃紺色の粒子は水粒子を表している．なお，見やすいように，表面にある水粒子は水色に着色した．粒子直径は 1.0mm であり，壁粒子 46,656 個，弾性体粒子 5,760 個，水粒子 27,648 個の合計 80,064 個の粒子を投入して計算を行った．弾性板の大きさは，x1 ，x2 ，x3 方向にそれぞれ 4mm，40mm，36mm の直方体型で，下壁との間に 2mm の隙間を設けた．弾性板の右側には，20×36 の断面に底面から 38mmの高さまで水が入っている．弾性体物性値は片持ち梁の場合と同じ，ヤング率 E = 1.859 MPa，ポアソン比 ν =0.49 とした．

Elastic plate

Water

ピークに達し，さらに元の形状へ戻りはじめている．梁の振動の場合と同様，ゲートの付け根と先端の応力分布
σ 22 の正負が逆になっており，ゲートがしなるように開閉していることがわかる．このような傾向は，初期水位や弾性板の長さといった条件が異なるものの，Carla ら 8)が実験及び計算によって観察した結果と一致している．
Figure 7 はゲートの自由端の初期位置から x1 ，x2 方向への変位を示している．ゲートは変形（反り）と復元を周期的に繰り返しており，水位の低下とともに次第にその振幅も小さくなっている．また，このときゲートに蓄えられたひずみエネルギーは，Fig. 8 に示す通り，ゲートが変形（反り）すると同時に蓄えられ，復元するとともに解放されている事がわかる．
結言
3 次元 SPH 法を用いて流体-弾性体の連成解析を行った．まず既存の弾性体解析モデルを 3 次元へ拡張し，片持ち梁の振動シミュレーションを行ったところ，理論解に近い結果が得られた．さらに，非圧縮性流体解析と組み合わせて連成解析モデルを構築し，弾性板のゲートを通過する液柱の崩壊現象の再現を試みた．その結果，ゲートの変形や液柱の崩壊の様子に関して，Carla ら 8) の計算及び実験と良く似た傾向のシミュレーション結果が得られた．
参考文献
1) 関根 英樹，飛行体の金属積層板に対する超高速衝突破壊現象の改良 SPH 法による数値解析，日本航空
宇宙学会論文集，Vol. 51(2003)，pp.45-51．

5.

2) 岡地 範明, SPH 法をベースにした非圧縮流れの計算手法の提案, 第 16 回流体力学シンポジウム, 講演要
旨集 (2002), pp.1-6.

Rigid wallFig. 2: Initial particle distribution of the dam-breakproblem with an elastic plate.
結果と考察
Figure 3 は時刻 t = 0s から 0.01s おきの計算結果を描画したものである．弾性体粒子の色は，カラーバーに示すように応力 σ 11 の大きさに比例している．なお，0.01
MPa 以下の応力は青，0.01 MPa 以上の応力は赤とした．
梁は周期 15.15Hz で 2 次元的に振動していることがわかる．また，梁が上向きに反ると，梁の付け根の上面側には圧縮力が下面側には引っ張り力が作用するので，梁の上面側には下向きの応力（青色）下面側には上向きの応力（赤色）が発生している．梁の振動にともなって，このパターンが交互に現れていることがわかる．梁の先端には付け根と逆のパターンが見られるが，これは梁がしなるように振動しているためである．このときの自由端の初期位置からの変位を Fig. 4 に示す．振幅の大きさは周期的に変化しながら次第に減衰していることがわかる．

4.

梁の長さ L を 40mm から 120mm まで変えた場合の振動周波数の変化を Fig. 5 に示す．図中の実線は，2 次元の基本モードに対する振動周波数の理論解である．長さ
L，厚さ H の片持ち梁の振動周波数は，√Hk 2Ef=(15)4π3ρ(1 − ν 2 )3) 渡辺 貴洋，非圧縮 SPH 法による二重拡散対流場の数値シミュレーション，日本機械学会 2006 年度年次大会講演論文集，Vol．2(2006)，pp.7-8.
4) 佐野 友哉，非圧縮性流体に拡張した SPH 法による自由表面および移動境界を有する流れの数値シミュ
レーション，日本機械学会流体工学部門講演会論文集，(2005)，pp.312.
5) 本郷 卓也，3 次元非圧縮 SPH 法における気液界面に作用する表面張力モデル，第 23 回数値流体力学シンポジウム，講演要旨集 (2009)，pp.165．
6) 伊藤 真澄，非圧縮ＳＰＨ法を用いたアーク溶融池内流れの数値シミュレーション ，日本機械学会年次大会講演論文集 ，Vol.2(2009)，pp.125-126
7) J. P. Gray, SPH elastic dynamics, Comput. Methods Appl. Mech. Engrg. , Vol. 190 (2001)，
pp.6641-6662．
8) Carla Antoci，Numerical simulation of ﬂuidstructure interaction by SPH, Computer and Structure，Vol．85 (2007)，pp.879-890．と与えられ，kL = 1.875 のとき基本モードとなる．両者を比べると，概ね近い値となっていることがわかる．

Figure 6 は弾性板のゲートを通過する液柱の崩壊現象の数値計算結果を初期ステップから 0.02s おきに描画したものである．ここで弾性体粒子の色付けは Fig. 3 と同様であり，σ 22 の値に対応している．はじめ静止していたゲートが流体によって押されて変形し，底面との間の隙間が広がるにつれて水が勢いよく流れ出していく様子がわかる．t = 0.06 ∼ t = 0.08s にかけてゲートの反りは
3c 2010 by JSFM
Copyright ⃝(a) t = 0.00s

x2 displacement of the free edge [mm]第 24 回数値流体力学シンポジウム
F10-4(b) t = 0.01s

Time [s]Fig. 4: Vertical position of the free edge of the cantilever.
(c) t = 0.02s(d) t = 0.03s(e) t = 0.04s(f) t = 0.05s(g) t = 0.06s

Frequency [Hz]: present computation(h) t = 0.07s

theoryσ 11[MPa]
Fig. 3: Oscillation of the cantilever from (a) t = 0.00 s
until (h) t = 0.07 s.

Length [mm]Fig. 5: Oscillation frequency of the cantilever.

4c 2010 by JSFM
Copyright ⃝(a) t = 0.000s

x1 and x2 displacement of the free edge [mm]第 24 回数値流体力学シンポジウム
F10-4(b) t = 0.020s

x1 displacement

x2 displacement

Time [s](c) t = 0.040s(d) t = 0.060s(e) t = 0.080s(f) t = 0.100s(g) t = 0.120s(h) t = 0.140sTotal strain energy [J]Fig. 7: x1 , x2 displacement of the free edge of the gate.

σ 22[MPa]
Fig. 6: Dam break problem with an elastic plate from(a) t = 0.000 s until (f) t = 0.140 s.

Time [s]Fig. 8: Total strain energy of the gate.

5c 2010 by JSFMCopyright ⃝