Fictitious Domain 法による着水衝撃の数値解析

Fictitious Domain 法による着水衝撃の数値解析

Numerical analysis of water-impact phenomena by the fictitious domain method

精密工学専攻

3

号 新谷 豪章

Hideaki Araya

1.

はじめに

着水衝撃とは物体が速度を持って水面に衝突するときに生じ る現象であり，船舶のスラミング，飛行艇の着水，宇宙船の海 上着水など例を見ることができる．荒海上を船舶が航行すると き，船首が波によって持ち上げられたのち海面にたたきつけら れることがある．これがスラミング

(slamming)

である．この ため船体外板などに局所的なダメージを被ったり，時には衝撃 圧が船首部甲板を座屈させるほどの大きさとなり海難事故を引 き起こすことさえある．このように衝撃圧を知るということは 船舶の安全性を確保するうえで必要である．このような物体の 着水現象をモデル化すると，固体と自由表面を有する液体の連 成問題になる．このような問題においては，液体領域を計算メッ シュで覆い，固体の運動や自由表面の変形に合わせてメッシュ を変形させたり，作り直しながら液体内の流れ計算を行うこと が多い．

ALE

^法

(Arbitrary Lagrangian-Eulerian method)

^に よる解析がその例である⁽¹⁾．しかし，このような方法では液 体領域に張ったメッシュを毎回作り直す作業が必要であり，計 算手順が複雑になるのが難点である．これに対して，液体領域 を含む広い領域に固定したメッシュを張り，その上で固体や自 由表面の位置を補足しながら流れ計算を行う方法がある．固体 と流体の連成問題に対する方法として

Fictitious Domain

法，

Ghost Fluid

法，

Immersed Boundary

法などが提案されてお り，自由表面問題に対する方法として，

Level Set

法や

VOF

法 などが提案されている．

室本⁽²⁾は，着水衝撃現象を固体と液体の連成問題としてで はなく，固体や液体の周囲に広がる気体領域も計算領域とする 気液固

3

相の連成問題として解析することを検討した．そして，

固相と気液相の連成解析には

Fictitious Domain

^{法を，気相と} 液相の連成解析には

Level Set

法を用いる計算手法を提案した．

さまざまな数値実験を通して，室本の方法は有効であることが 示されたが，実用解法とするのには計算精度の向上を図る必要 があることがわかった．そこで，本研究では，室本の計算方法 の精度を向上させることを目的として行う．

2.

流れの支配方程式

Fig.

１のような

2

次元のモデルを考える．気相と液相内の流 体は非圧縮性粘性流体とする。このとき流れの支配方程式は次 の非圧縮ナビエ・ストークス方程式と連続の方程式である．

∂u

∂t + (u · ∇)u = 1 ρ ∇ · ff + F in Ω

^′

(1)

∇ · u = 0 in Ω

^′

(2)

である．ここで，

Ω

^′は気相と液相を合わせた流体領域を表す．

また

u

^は速度，

ff

^{は応力テンソル，}

ρ

^は密度，

F

^{は外力を表す．}

gravity

liquid gas

solid

Fig. 1: A computational model for the water entry ploblem of a rigid body

外力として重力を考える．流体をニュートン流体として，

ff

を 次のように与える．

ff = −pI + µ [

∇u + (∇u

^T

) ]

(3)

ここで，

p

^は圧力

, µ

は粘性係数を表す．固体表面を

∂ω

^する．

そのとき，固体表面境界

∂ω

^上に

u = u

b

on ∂ω (4)

を課す．ここで

u

bは固体の速度である．

Fig.1

の

4

辺の壁境 界には滑り条件を課す．

3. Fictitious Domain

法

Fictitious Domain

法では，流体が占める領域

Ω

^{′を計算領域} とするのではなく，

Ω

^{′に固体が占める領域}

ω

^を加えた

Ω(Fig.1

の長方形領域の内部全体

)

を計算領域とする．そして，固体境界

∂ω

上の境界条件を，境界条件としてではなく，

Ω

^内の

u

^に対 する拘束条件として用いる．式

(1), (2)

^に

Ficititious Domain

法を適用して弱形式を導くと，次式のようになる．

ρ

∫

Ω

[ _∂u

∂t + (u · ∇)u ]

· u

^∗

dΩ −

∫

Ω

p∇ · u

^∗

dΩ + µ

∫

Ω

[ ∇u : u

^∗

+ ∇u : (∇u

^∗

)

^T

] dΩ

− ρ

∫

Ω

F · u

^∗

dΩ =

∫

ω

– · u

^∗

dω (5)

∫

Ω

p

^∗

∇ · u dΩ = 0 (6)

∫

ω

–

^∗

· (u − u

b

) dω = 0 (7)

ここで，

–

はラグランジュ未定乗数であり，アスタリスク

∗

のついたものは重み関数を表す．

式

(7)

^は式

(4)

^{から導かれる式}

(5)

に対する付帯条件である．

このように，ナビエストークス方程式の弱形式に

–

^{の境界積分} 項を付加することにより，物体表面境界

∂ω

^{上での境界条件を} 拘束条件化する事ができ，物体の存在を意識せずに

Ω

^全体を 計算領域と考えることができる．そのため，物体の移動に伴う メッシュの再構築の必要がなくなる．

また

Fictitius Domain

法ではラグランジュ未定乗数

–

^を

ω

に沿って積分することで流体力が簡単に求まるため

,

流体構造 練成問題には非常に有用性が高い．

4.

弱形式の離散化

4.1

空間方向の離散化

空間方向の離散化には有限要素法を用いる．有限要素として 三角形要素を用いる．速度は

Fig.2(a)

のように，要素を４つの 小三角形に分割し，各小三角形内で

x, y

の１次式で近似する．

圧力は

Fig.2(b)

のように全体で

x, y

の

1

次式で近似する．

領域式

ω

を覆う有限要素メッシュとは別に，

(5),(7)

の

ω

内 の積分のために

ω

を覆うメッシュを別に用意する．要素は三角 形とし，

–,–

^˜は要素内で

x, y

の

1

次式で近似する．以上より，

式

(5)-(7)

^{は離散化され，}

M dU

dt + [A(U) + D] U − HP − F − BM

ω

Λ = 0 (8)

H

^T

U = 0 (9)

B

^T

U = U

b

(10)

となる．ここで，

M

時間微分項，

A(U)

は移流項，

D

は拡散項，

H

は圧力項，

BM

ωは

–

の項に対応する係数行列である．

M

の上付きバーは対角化を意味する．また，

F

は外力項を表す．

(a) velocity (b) pressure Fig. 2: Nodes of a Bercovier-Pironneau element

4.2

時間方向の離散化

時間方向の離散化には差分法を用いる．時間軸を一定の長さ

∆t

の小区間に分割し，時刻

t

ⁿ

= n∆t

^と時刻

t

ⁿ⁺¹

= (n +1)∆t

に挟まれた区間を考える．そこで，式

(8)-(10)

^{を時間に関して} 離散化するとともに，次のような

3

組の方程式群に分解する．

M U

^∗

− U

ⁿ

∆t + [

A(U

ⁿ

) + D ]

U

ⁿ

− HP

ⁿ

−F − BM

ω

Λ

ⁿ

= 0 (11)

 



 

M U

^∗∗

− U

ⁿ

∆t − H (

P

ⁿ⁺¹

− P

ⁿ

)

= 0 H

^T

U

^∗∗

= 0

(12)

 



 

M U

ⁿ⁺¹

− U

^∗∗

∆t − BM

ω

(Λ

ⁿ⁺¹

− Λ

ⁿ

) = 0 B

^T

U

ⁿ⁺¹

= U

b

(13)

ここに，上付き添字

∗, ∗∗

は中間的な速度を意味する．時刻

t

ⁿ における物理量に値を知って，時刻

t

^{n+1における値を求める} ために，式

(11)-(13)

を次のステップで段階的に解いていく．

step1:

^中間流速

U

^∗を式

(11)

^{より求める．}

step2:

^式

(12)

^の第

1

^式を

U

^∗∗について解き，その結果を式

(12)

^の第

2

^{式に代入すると，}

( H

^T

M

⁻¹

H )

Φ = − H

^T

∆t U

^∗

(14)

を得る．ここに，

Φ = P

ⁿ⁺¹

− P

ⁿ

(15)

である．式

(14)

を解いて圧力修正量

Φ

を求め，式

(15)

よ り

P

ⁿ⁺¹を計算する．そして，式

(12)

の第

1

式によって

U

^∗∗を求める．

step3:

式

(13)

の第

1

式を

U

ⁿ⁺¹について解き，その結果を 式

(13)

の第

2

式に代入すると，

∆t(B

^T

M

⁻¹

BM

ω

)Ψ = U

ⁿ⁺¹_b

− B

^T

U

^∗∗

(16)

を得る．ここに

Ψ = Λ

ⁿ⁺¹

− Λ

ⁿ⁺¹

(17)

である．式

(16)

^を解いて

Λ

^の修正量

Ψ

^{を求め，式}

(17)

より

Λ

ⁿ⁺¹ を計算する．そして，式

(13)

^の第

1

^式によっ て

U

^{n+1を求める．}

連立

1

^{次代数方程式}

(14),(16)

の求解には，共役勾配法を用いる．

5.

レベルセット法

自由表面の位置を特定するために

,

レベルセット法を用いる．

流体領域

Ω

^{′において関数}

φ(x, y, t)

^{を定義する．}

φ

^は

φ

{ > 0 (

液相において

)

< 0 (

気相において

) (18)

とし，

|φ|

は自由表面 からの距離を表すとする．このとき，

φ = 0

の等値線が自由表面を表すことになる．そこで，各時刻におい て

Ω

^′内の

φ

の分布の中から

φ = 0

の曲線を見つければ自由表 面の位置と形状を把握することができる．

φ

の時間変化は移流 方程式

∂φ

∂t + (u · ∇)φ = 0 (19)

に従う．式

(19)

^を解いて

φ

を求め続けていくと，計算誤差の影 響によって距離関数を表すという性質が失われてしまう．

φ

^が 距離関数であるという性質を保つための処理が行われる．再初 期化は

∂φ

∂τ + S(φ)(|φ| − 1) = 0 (20)

を解くことで行われる．

τ

は仮想時間である．また

S(φ)

^は

S(φ) = √ φ

φ

²

+ (9h)

²

(21)

で定義される．流体内の各点での密度と粘性係数は次式で計算 する．

ρ(φ) = (1 − H (φ))ρ

l

+ ρ

g

· H (φ) (22) µ(φ) = (1 − H(φ))µ

l

+ µ

g

· H(φ) (23)

ここで，

ρ

l

, ρ

gはそれぞれ液体，気体での密度，

µ

l

, µ

g はそれ ぞれ液体，気体の粘性係数である．

H(φ)

は次式で定義される 関数である．

H(φ) =

 

 

 

 

 



1 : φ > α

0 : φ < −α

φ + ϵ 2ϵ

1 2π sin( πφ

α ) : |φ| ≤ α

(24)

ここで，

α

^{は要素の代表長さを}

h

^{としたとき，}

α = 2h

^の大 きさである．式

(24)

^{は自由表面を厚さ}

4h

^{の幅をもって帯状の} 領域で表すことを意味している．

6.

計算手順

6.1 3

相計算のアルゴリズム

3

相計算のアルゴリズムは以下のとおりである．

1.

^式

(22)

^，

(23)

^{，節点ごとに密度}

,

粘性係数の値を計算する．

2.

流れ計算を行い，速度

u,

圧力

p

ラグランジュ未定乗数

–

を求める．

3.

^{移流方程式}

(19)

を解き，レベルセット関数

φ

^{の新たな分} 布を決定する．

4.

式

(20)

を解いて

φ

の再初期化を行う． 

5.

固体に作用する流体力を計算する．

6.

求めた流体力を基にして固体の運動方程式を解き，固体の 速度と変位を求める．

7.

固体の変位に合わせて固体領域

ω

^{のメッシュの位置を更} 新する．

8.

時間を

∆t

だけ進めて手順

1.

から繰り返す．

6.2

流体力の計算

Fictitious Domain

法において，ラグランジュ未定乗数

–

は

– = ff · n (25)

という意味を持つ．ここに

n

は固体表面に立てた外向き単位法 線ベクトルである．そこで，

F

λ

= −

∫

ω

–dω (26)

によって，固体が流体から受ける力を計算することができる．し かし，

F

λには浮力の効果が反映されないので，浮力は別に計 算しなければならない．そこで，浮力

F

bを

F

b

= −g

∫

ω

ρdω (27)

で計算する．以上より，固体に作用する流体力

F

^は

F = F

λ

+F

b

で与えられる．

6.3

固体の運動方程式

固体の運動方程式を次式で与える．

M d

²

x

dt

²

= M g + F (28)

ここに，

x

^{は物体の重心の座標，}

M

^{は物体の質量，}

g

^は重力等 の外力，

F

^{は流体力を表す．式}

(28)

^{の解法には}

4

^{次のルンゲ・}

クッタ法を用いる．

7.

メッシュの細分化

室本は

Cheng

ら⁽³⁾が行った着水実験と比較し，計算精度の

検証を行った．そのモデルは

Fig.4

に示す．しかし，時々刻々変 化する物体の位置や物体に加わる圧力を正確に捉えるにはまだ 改良が必要であった．

Loon

^ら(4)は，物体周辺の有限要素メッ シュを細分化することによって計算精度が向上することを報告 している．しかし，彼らは，細分化の際にコントロールポイント が節点に一致するように物体形状に合わせてメッシュを生成して いた．物体形状に合わせてメッシュを分割するということはラグ ランジュ的解法である

ALE(Arbitrary Lagrangian-Eulerian)

法と同じであり，物体の形状や位置とは独立にメッシュを張る ことができるという

Fictitious Domain

^{法の長所を損なうもの} であると考えられる．よって，時々刻々物体が移動する本問題 に対して有効な手段であるとは言えない．そこで，本研究では 物体近傍のメッシュに対して細分化を施すことにより計算精度 が向上するのではないかと考えた．

Fig.3

を用いて細分化の手 順を説明する．細分化前のメッシュは

Fig.3(a)

のように規則的 に要素が並んだ構造とする．はじめに

Fig.3(b)

の太枠に囲まれ た領域のように物体が存在する領域より大きく細分化領域をと る．そして，

Fig.3(c)

のように細分化領域内の三角形要素を

4

分割する．このままでは細分化領域とそうでない領域の

2

つの 領域に分かれてしまう．最後に

Fig.3(d)

^{に示すように，細分化} 領域の境界上にある外側の三角形要素を

2

^{分割する．また，時々} 刻々物体の移動に合わせて細分化を行うのではなく，物体があ る一定の距離を移動した場合，そのときの物体の位置に合わせ，

Fig.3(a)

の初期のメッシュの状態から改めて細分化を行う．

mesh

solid

(a) Original mesh (b) Area of subdivision

(c) Division into four (d) Division into two Fig. 3: The way of subdivision of mesh

8.

矩形物体の着水解析

Cheng

らが行った矩形物体の着水実験と比較し，計算精度

の検証を行う．矩形物体の幅

L

を

0.076m

とし，鉛直下向きに

0.33m/s

の初速度を与える．また，矩形物体の質量を

10.22kg

とする．液体の密度と粘性を

10

³

kg/m

³

,10

⁻³

Pa·s

とし，気体の 密度と粘性を

10kg/m

³

,10

⁻⁴

Pa·s

とする．今回用いるメッシュ の要素数は

23232,

^節点数は

11771

である．細分化用に用いた物 体領域

ω

^{メッシュは要素が}

1200,

^節点が

641

^{である．また，細} 分化を行わない場合に用いた固体領域

ω

^{メッシュは要素が}

432,

節点が

241

^である

. Fig.5

に静止液面から測った矩形物体の落下 距離の時間変化を示す．

Cheng

らの実験値とはやや離れてはい るがメッシュの細分化を行い計算を行った方がより実験値に近 づくことが確認できる．

Fig.6

は矩形物体底面の中心に働く圧 力の時間変化を示している．メッシュの細分化を施すことによ り，圧力の時間方向における圧力の振動を抑えられていること がわかる．また，圧力の値も向上し，やや実験値に近づいてい る．しかし，水面に着水瞬間におけるピーク圧が細分化を行っ た場合は

0.89kPa,

行わなかった場合は

0.68kPa

という実験値 とは大きくかけ離れた結果となってしまっている．

gravity

20L L

10L 15L 3.14L

Fig. 4: Water entry of a rectangular body

9.

おわりに

固体の着水現象を気液固

3

相連成問題として解析する計算手 法の精度向上を目指して，固体の移動に合わせて固体周辺の要 素を細分化する機能を付加し，その効果を検証した．実験結果 との比較を通じて，要素細分化の有効性を確認した．しかし，

着水の瞬間に作用する衝撃圧の大きさについて，計算値と実験 値の間には

5kPa

ほどの大きな差があり，今後のさらなる検討 が必要である．

Time[s]

Displacement[cm]

Experiment(Cheng  et  al.,  1974) Present  caluculation(with  subdivision) Present   caluculation(without  subdivision)

Fig. 5: Time history of the displacement of the body

Time[s]

Pressure[kPa]

Experiment(Cheng   et   al.,   1974)

Present   caluculation(without  subdivision) Present  caluculation(with  subdivision)

Fig. 6: Time history of the pressure acting at the center of the impact surface

参考文献

(1)

^{田中 聖三}

,

自由表面流れ問題に対するバックグラウンドメッ シュを用いた

ALE

有限要素法の構築研究， 博士論文， 中 央大学

(2006).

(2)

室本 悠介

, Fictitious Domain

法とレベルセット法の併用 による多相連成解析， 修士論文， 中央大学

(2008).

(3) Cheng, R. Y. K. and Leland, T. J. W. , Numerical so- lution for low-velocity penetration of rigid body into still fluid ”, Numerical Methods in Fluid Dynamics (edited by C. A. Brebbia and J. J. Connor), Pentech Press, London, 1974, pp. 272-289

(4) R. van Loon, P. D. Anderson, A combined fictitious domain/adaptive meshing method for fluid-structure in- teraction in heart valves”, Int. J. Numer. Meth. Fluid, 46 (2004), pp.533-544

(5) R. Glowinski, T. W. Pan, T. I. Hesla, D. D. Joseph and J. Periaux, A Fictitious Domain Approach to the Direct Numerical Simulation of Incompressible Viscous Flow past Moving Rigid Bodies: Application to Particu- late Flow”, J. Comput. Phys. , 169(2001), pp. 363-426 (6) M. Sussman, P. Smereka and S. Osher, A Level Set Ap-

proach for Computing Solutions to Incompressible Two-

Phase Flow”, J. Comput. Phys. ,114(1994), pp. 146-159