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

非数値拡散フィルタの気液混相流数値計算への適用性の検討
Application of Anti Numerical Diﬀusion Filter to Gas-Liquid Multiphase Fluid Simulation
藤岡奨, 水資源機構, 埼玉県さいたま市中央区新都心 11-2, susumu fujioka@water.go.jp

牛島省, 京都大学学術情報メディアセンター, 京都府京都市左京区吉田本町, ushijima@media.kyoto-u.ac.jp

Susumu Fujioka, Japan Water Agency, 11-2 Shintoshin,Chuo-ku,Saitama, Japan
Satoru USHIJIMA, Kyoto University, Yoshida-honmachi, Sakyo-ku, Kyoto 606-8501, Japan
In volume of ﬂuid (VOF) method, the accuracy of advection of VOF function is essential to predict highly deformedinterface between the liquid and gas. An anti numrical diﬀusion (AND) ﬁlter has been proposed to overcome
numerical diﬀusion problem in advection. The accuracy of AND ﬁlter has been conﬁrmed through various test
cases. In this paper, the AND ﬁlter is applied to 5th order TVD scheme in advection of VOF function. Theeﬀectiveness of AND ﬁlter in multi-phase ﬂuid simulation is conﬁrmed through test cases of sloshing. The results
showed that AND ﬁlter is eﬀective to keep sharp interface between the ﬂuids.

緒言
気液混相流れの数値形跡手法の一つである VOF 法 1においては，VOF 関数値によって界面を識別するため，幾何学形状にて界面形状を表現する手法に比べて簡便に複雑な水面形状に追随することができる．一方で，VOF関数の移流計算を精度よく解かなければ数値拡散によって界面幅が広がり，計算の進行とともに界面形状の予測精度が落ちる．この問題に対して，種々の対策が研究されている 2, 3, 4 ．
既往文献 5, 6 においては，VOF 関数の移流計算の精度向上をさせるために非数値拡散フィルタ (Anti NumericalDiﬀusion Filter, AND フィルタ) が提案されている．ANDフィルタは，STAA 法 7, 8 と同様に移流方程式を解いた後のスカラー分布へ適用し，数値拡散を補正する働きをするものである．
AND フィルタ導入の利点としては，移流計算に伴うVOF 関数の数値拡散を抑えて急勾配を保つことができることに加え，スカラーの微分値等を使用しないためコンパクトであり，アルゴリズム上２次元から３次元への拡張が容易であること，並列計算への適用が比較的容易であることが挙げられる．更に，隣接セル間での値のやり取りを基本とするアルゴリズムであるため，非構造格子への適用も期待できる．文献 9 においては，並列化した流体計算への AND フィルタの適用が検討されているが，界面厚さ保持の性能評価については目視による評価にとどまっている．本報告では，並列化された流体計算に AND フィルタを適用し，理論解と比較して混相流れにおける界面厚さの定量的評価を行い，AND フィルタの有効性を確かめた．
混相流流体計算法として MICS(Multiphase Incompressible ﬂow solver with Collocated grid System ) 10, 11を用いた．VOF 関数の移流スキームには，既往研究 6 にて AND フィルタとの組み合わせにより高精度の移流計算が可能であるとされた高次精度の TVD スキーム 12, 13を用いた．

1.

2. 計算手法(1) 混相流れ数値計算法
本研究では，MICS(Multiphase Incompressible ﬂow
solver with Collocated grid System ) 10, 11 を用いて混相流流体計算を行った．基礎式は以下に示す質量保存則，非圧縮性条件式，及び運動方程式である．
∂ρ∂+(ρuj ) = 0
∂t
∂xj
∂uj=0
∂xj

(1)
(2)

1

∂
1 ∂p
∂ui+(ui uj ) = fi −∂t∂xj
ρ ∂xi]∂1 ∂ [ ∂(µui ) +(µuj ) · · · (i = 1, 2, 3)+ρ ∂xj ∂xj
∂xi(3)

ここに，t は時間，xi は３次元直交座標系の座標成分である．ρ, µ, p は計算セルにおいて体積平均された密度，粘性率，圧力である．ui は計算セルにおいて質量平均された xi 方向の流速成分である．fi は xi 方向の外力である．式 3 をコロケート格子上で有限体積法により離散化し，MAC 系解法に基づいて予測段階，圧力計算段階，修正段階を経て次時間ステップの流速及び圧力を更新する．
予測段階においては，C-ISMAC 法による完全陰解法を用い，圧力計算段階においては，C-HSMAC 法を用いた．次時間ステップにおける流速及び圧力を更新した後，VOF関数の移流計算を５次精度の TVD スキームあるいは 1次風上差分法を用いて行った．上記移流計算後，非数値拡散フィルタを適用し，次時間ステップにおける VOF 関数を得た．
(2) 非数値拡散フィルタの概要
Fig.1 のフローチャートに示される様に，AND フィル
タは VOF 法におけるスカラー関数の移流計算が終わった段階で適用され，移流計算に伴う数値拡散を補正する働きをする．具体的には，周辺の計算セルよりも値の小さなスカラー値をもつセルからスカラー最急勾配方向に修正量を移動させ，移動操作後のスカラー値が 0 と 1 の間に収まる範囲内でスカラー勾配を急峻化させるアルゴリズムである．各段階の手順の詳細は，文献 6 にて検討されている．以下では，計算手順を概説する．
移流計算を終えた VOF 関数 φ に対して前処理を適用した後，修正量の移動方向を Fig.2 に示される様に直線上に連続して隣接する４点を用いて定める．移動ベクトルの定義に必要な 4 点は以下の通りである．

Z 点 スカラー値が 0 の点である．
F 点 Z 点に隣接しており，スカラー値が正である点である．
1,2 点 Z 点から F 点の方向に延長した移動ベクトルの先に位置する点である．
隣接する 4 点の始点は VOF 関数が 0 となる点 (Z 点)であり，始点を除く 3 点の VOF 関数が全て正の値を持っていた場合，中間の２点に移動ベクトル s̃(a) を定義する．
ただし，a は s̃ の始点となるセル中心の空間座標を示す配列である．また，s̃ の終点側となるセル中心の空間座標を示す配列を b とする．
これら４点から移動ベクトル s̃(a) を定義した後，Fig.2に示される様に，F 点を挟んで Z 点と 1 点の中央差分にCopyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-5s(a)VOF 関数の移流計算

C

非数値拡散フィルタの適用
1

前処理

2変数の初期化

3

移動ベクトルの定義

4

スカラー移動量の定義

5スカラー移動量補正

6スカラー移動実行

B
AΔφ0.0Z

Fx=a2

1x=bFig. 2: 移動ベクトル s(a) の定義

ただし，T は式 10 で定義される．

T = {j ∈ N |j ≤ M (b)}Fig. 1: AND フィルタ計算手順

ここに，φ(b) は移動ベクトル s の終点側セルの VOF 関数値であり，当該セルへ M (b) 個の方向から移動方向ベクトル s が定義されているものとする．したがって，修正量 r(a) は式 11 によって定義される．てスカラー勾配 ∆φi (a) を式 4 にて算出する．

∆φi (a) =φ(a + s̃i (a))
2|s̃i (a)|rk (a) = min(r̃(a), rmax (a + s(a)))(4)∆φj (a) = max ∆φi (a)i∈S(5)φ (a) = φ(a) −i=1

(6)

3.

ここに，N は自然数を表す集合であり，Ñ (a) は a を始点として定義される s̃i (a) の総数である．また，sk (a) に沿って移動する修正量の初期値 r̃k (a) を式 7 によって定義する．

r̃k (a) =φ(a)
N (a)

(7)

j∈T

ri (a) +rj (a)

(12)

j=1

数値計算結果∆η(t, x) =Hcos(kx) cos(ωt)2(13)

ここに，H は波高，k は波数，ω は角周波数，t は時刻，x は波の進行方向に沿った位置座標である．波長 L に対して水深 h が十分に大きい場合深水波呼ばれ，以下の近似が成り立つ．
√gL(14)c=2π(8)

従って，式 7 によって初期化された修正量 r̃k (a) を式 9によって定義される rmax によって制限し，修正操作後のセルの VOF 関数が 1 を超えないように調整する．

rmax (b) = min(min(r̃j ), (1 − φ(b))/M (b))∑M (a)

3.1 完全流体におけるスロッシング
AND フィルタの気液混相流れへの適用性を確かめるために，矩形水槽内の完全流体のスロッシングを数値計算し，理論解と比較することによって適用性を確かめた．
VOF 関数の移流計算には 5 次の TVD スキーム及び１次風上差分を用い，それぞれのケースに対して AND フィルタを適用し，その効果を確かめた．微小振幅波理論によると，波高が水深に比べて十分に小さく流体の粘性が無視できる場合，静水位からの水位変位 ∆η は以下の式で表現される．

N (a) は Ñ (a) と同様に a を始点として定義される s(a)の総数である．
修正量 rk (a) の算出には始点 a における VOF 関数値
φ(a) および，移動先のセルの VOF 関数値 φ(a + s(a)) 及
び他のセル a∗ から同じセルへ移動予定の修正量 r(a∗ ) の値を考慮する必要がある．ただし，a，a∗ の関係は式 8にて表される．
b = a + s(a) = a∗ + s(a∗ )∑N (a)∗

ただし，集合 S は式 6 で定義する．

S = {i ∈ N |i ≤ Ñ (a)}(11)

このように修正量を調整する操作によって空間的にバランスのとれた修正操作を行うことが可能である．また，以下では移動元セル a における rk (a) を移動先セル b で参照する場合には，rm (b) と表現するものとした．
以上の手順により算出された s(a) および r(a) を用いて，フィルタ操作完了後のスカラー φ∗ を式 12 にて求める．

F 点 a を始点とする移動ベクトル s̃i (a) のうちで，スカラー勾配 ∆φi (a) が最大値と等しくなる方向を最終的に移動ベクトル sk (a) として選出する．このため，a に対して複数の sk (a) が定義されることもある．式 5 に s(a)の定義式を示す．
sk (a) = s̃j (a) where(10)(9)ここに，g は重力加速度であり，c は波速である．

2

Copyright c 2010 by JSFM

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

以下では一辺 0.5×0.5×1.0[m] の直方体水槽内に生じる微小振幅波を対象として検証計算を行った．以下の検証計算では上記の水槽を 16 × 16 × 32 のセル分割として計算した．
初期状態 t = 0 では水槽に静水位 h =0.5[m] の液体が入っているものとしたため，位置 x における水位 η(0, x)は次式で表される．

Hcos(kx)η(0, x) = 0.5 +2は界面厚さが計算セル幅と比較してどの程度まで広がっているかを示しており，界面幅がセル幅より小さくなるケースにおいては最小値 0.9 となる．

w(b) =(18)

ここで，zα (b) は xy 平面上において，b = (i, j) の位置の計算セル線形内挿によって VOF 関数が α の値を持つ高さである．
また，式 18 にて定義される界面厚さ w(b) の平均値wave. を算出し，計算領域全体平均界面厚さを比較することによって AND フィルタの効果を確かめた．

(15)

ただし，波高 H = 0.1m, 波長 L =1m，波数 k = 2π とした．
流 体 の 粘 性 を 全 て ０ と し ，液 体 の 密 度 を 1.0 ×103 kg/m3 ，気体の密度を 1kg/m3 とした．

∑zwave. =

500

500

z0.05 (b) − z0.95 (b)∆zS w(b)Nx Ny(19)∑
ここに， S は水平面全体で和を取ることを示し，Nx , Ny
はそれぞれ x, y 方向のセル数である．
VOF
1.0
VOFm
0.95
VOFm+11000

1.00.0z0.95mhyxdz

m+1m+n-1z0.05m+nFig. 4: Deﬁnition of interface width以上の条件にて t = 0 から t = 4[s] まで数値計算を実行し，得られた誤差を Fig.5 に示す．
AND フィルタを使用したケースに比べて，AND フィルタを使用しないケースは時間経過とともに誤差が増加している．これは，数値拡散による影響で界面幅が広がったためと考えられる．一方で，AND フィルタを使用したケースでもある程度の誤差が発生しているが，誤差の上限が抑えられている．ただし，AND フィルタを使用した
場合，誤差の増加パターンが単調ではなく，周期的に増減を繰り返しながら緩やかに増加している傾向が特徴的である．
界面厚さについては Fig.6 に示されるように，ANDフィルタを適用したケースにおいては得られた界面厚さが薄く保たれていることが確認できる．
x = 0 において界面厚さの上限と下限となる z0.05 ，z0.95及び界面の中心線となる z0.5 の位置の時間変化を Fig.7に示す．5 次 TVD スキームに AND フィルタを適用したケースでは，AND フィルタの効果によって界面の拡散が防がれていることが確認できる．1 次風上差分に ANDフィルタ適用したケースでは，界面幅を薄く保つ効果とともに，位相のずれを修正していることが確認できる．

GravityFig. 3: Outline of Rectangular Tank
誤差を評価するため，時刻 t における VOF 関数の理論解 VOFe (a) を式 16 にて定義した．
top
 1.0 · · · η(t, x) ≥ z (a)VOFe (a) =0.0 · · · η(t, x) ≤ z bot (a)(η(t, x) − z bot (a))/∆z · · · otherw

VOFm+n-10.05
VOFm+n0.0 z(16)ここに，a はそれぞれ x, y, z 方向のセル番号を格納した配列であり，各軸方向のセル番号をそれぞれ ix , iy , iz とした場合，a = (ix , iy , iz ) である．z top (a), z bot (a) はそれぞれセル a の z 軸方向正方向，負方向のセル境界の z 座標値である．∆z は，z top (a) と z bot (a) の差である．式 16にて定義された VOFe (a) を用いて誤差 Error を式 17 にて定義した．
∑Error =|VOFe (a) − VOF(a)|∆V (a)(17)V∑
ここに， V は計算領域全体で和を取ることを意味す
る．VOFe (a) は，当該セルの理論的 VOF 関数値であり，VOF(a) は，数値計算結果から得られた VOF 関数値である．∆V (a) は，当該セルの体積である．
また，界面厚さ w(b) は Fig.4 に示される様に計算セルの VOF 関数値の線形空間内挿を元に，VOF 関数の 0.05
及び 0.95 の等高線の間隔にて定義した．界面厚さ w(b)33.2 粘性流体におけるスロッシング
次に，外力によって周期的な変位を生じる３次元矩形水槽のスロッシングの数値計算を行い，実験値 14 との比較を行った．水槽の緒元は幅 1.0[m]，奥行 0.2[m]，高さ 1.0[m]であり，静水位が 0.5[m] である．振幅 A =0.0093[m]，角速度 ω =5.311[rad/s] の条件で水槽を振動させ，水槽内に生じるスロッシング現象を数値解析対象とした．外力
加速度 f は式 20 で与えられる．
f = Aω 2 sin ωt(20)液相の密度を 1.0 × 103 [kg/m2 ]，動粘性係数を 1.0 ×
10−6 [m2 /s] とし，気相の密度を 1.0[kg/m2 ]，動粘性係

Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-5*+,(*+,(-./,%"#!"!!!&%!"!!!%$"#./012

01121

!"!!'

!"!!!$!"!!!#!

()*#()*#+,-*$!"#!

!"(''"(##"())"(!

$&&"#

!

*345678(a) 5th Order TVD scheme!"# $ $"#
(342567%%"#'(a) 5th Order TVD scheme

!"!!'

#"
()*#+,-.
()*#

!"!!&

()*#+,-.
()*#

*/012

#!

/0010

!"!!%

!"!!$
"
!"!!#

!

!
!

!"'

##"' $ $"'2345678

%%"'

&!

!$"

##$"%%$"&&$"'3452678(b) 1st Order up wind(b) 1st Order up wind

Fig. 5: Time histories of Errors (inviscid ﬂuid)

Fig. 6: Time histories of mean width of interface (inviscid ﬂuid)

数を 1.0 × 10−5 [m2 /s] とした．格子分割は 40 × 8 × 40 とした．Fig.8 に実験値及び x = 1.0 における水位の時系列を示す．Fig.8(a) は VOF 関数の移流計算に 5 次 TVDスキームを使用した計算例であり，Fig.8(b) は VOF 関数の移流計算に 5 次 TVD スキームを使用した後，ANDフィルタを適用した計算例である．AND フィルタを使用しないケースでは実験値との良好な一致を見せているが，AND フィルタを使用したケースでは 4.0[s] の時点までは実験値と良好な一致を見るが，4.0[s] 以降最高水位の値が過小評価されている傾向が確認された．
Fig.9 に界面厚さの計算領域平均 wave. の時系列を示す．
Fig.9 によると AND フィルタを適用したケースでは，界面幅が相対的に薄く保たれていることが確認できる．

波数 ω はそれぞれ 0.66π[s−1 ], 1.32π[s−1 ] とした．

f = Aω 2 sin ωt(21)

時間刻み ∆t を 0.01[s] とし，10[s] まで数値計算を行った．

Fig.10 に数値計算結果を示す．複雑な水平方向加速度によって大きく変形した界面を安定して追跡できていることが確認された．
4.

結言

最後に，外力によって複雑な変位を生じる３次元矩形水槽のスロッシングの数値計算を行なった．水槽の緒元
は 0.3 × 0.3 × 0.14[m] とし，計算領域を 60 × 60 × 28 のセルに分割した．流体計算においては計算領域を 8 分割し，並列計算を行った．並列計算における計算領域境界部のりしろ部分として 4 列分の格子を用いた 9 ．液相及び気相の物性は検証計算と同様とした．水槽に与える外力加速度は式 21 の形によって与え，静止時の水位を 0.1[m]とした．x 方向，y 方向の振幅 A は 1.5[m],0.09375[m], 周

本研究においては，AND フィルタの気液混相流れ数値計算への適用性を確認するために，自由水面変動を伴う完全流体及び粘性流体の数値計算を行ない，AND フィルタによって界面厚さを保持できることを確かめるとともに，理論解及び実験値との比較によって手法の妥当性を確認した．その結果，AND フィルタの効果によって，界面厚さが一定に保たれることが示された．一方で，粘性流体のスロッシング再現計算においては，部分的に界面位置の予測精度が落ちる例も確認された．界面位置の予測精度の向上については，原因の調査も含め今後の検討課題としたい．さらに，３次元矩形水槽内で大きな液面変形を伴うスロッシング現象に本手法を適用し，水槽内の複雑な液面形状を安定して数値計算できることを確認した．

4

Copyright c 2010 by JSFM

3.3

3 次元矩形水槽スロッシング第 24 回数値流体力学シンポジウム
E10-5

t=0.1 [s]t=0.2 [s]t=0.3 [s]t=0.4 [s]t=0.5 [s]t=0.6 [s]t=0.65 [s]t=0.7 [s]Fig. 10: Sloshing in a rectangular parallelepiped tank5

Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-5)*+$,*-./!"!$,)*+$,*-./!"$,)*+$,*-./!"0$,)*+$123+,*-./!"!$,)*+$123+,*-./!"$,)*+$123+,*-./!"0$,!"%$!")>8;?@/A?B?A/69=4,56/!7!"$$!"$!"#$!"#,-.%/-012!"!%,-.%/-012!"%,-.%/-012!")%345"/60789:;:</+))*=!"(!"%!"'!"&!"%!"$!

!"$&&"$''"$(("$!"##!

*)89:5;7(+!

(a) 5th Order TVD scheme+,-)./012!"!%234!
+,-)./012!"%234!
+,-)./012!"5%234!
+,-)2!"!%234!
+,-)2!"%234!
+,-)2!"5%234!

!"),-.%/01.2-345!"!%,-.%/01.2-345!"%,-.%/01.2-345!")%678"293:;<=>=?2+))*@!"(A;>BC2DBEBD29<@!"'6!"&!"%!"$!"#&,C9?6D=(a) 5th Order TVD scheme!"($!"'!"&!"%!"$!

!"%))"% * *"%789:;<=##"%!"#$!

*$&(+!

,F<B9G@(b) 1st Order up wind(b) 1st Order up wind

Fig. 7: Time histories of water level (inviscid ﬂuid)

Fig. 8: Time histories of water level (viscous ﬂuid)

参考文献

schemes based on rational functions. J. GeophysicalResearch-Aomosphere, Vol. 107, No. D22, 11 2002.

1. C. W. Hirt and B. D. Nichols.

Volume of ﬂuid

method for the dynamics of free boundaries.
Comp. Phys., Vol. 39, pp. 201–225, 1981.
2. Yokoi K.

Art. No. 4609.

J.

5. 藤岡奨, 牛島省. 気液混相流計算のためのスカラー交換型非拡散フィルタの開発. 応用力学論文集, Vol. 13,Eﬃcient implementation of THINC

scheme: A simple and practical smoothed VOF algorithm. J. Comput. Phys., Vol. 226, pp. 1985–2002,
2007.

pp. 271–280, 2010.
6. 藤岡奨, 牛島省. VOF 関数の移流計算における数値拡散抑制手法の検討. 計算力学論文集 (掲載予定).

3. Honma Y. Xiao F. and Kono T. A simple algebraic
7. 池端昭夫, 正肖鋒. 保存型自由界面捕獲スキームと固気液三相流への適用. 日本機会学会年次大会講演論文

interface capturing scheme using hyperbolic tangentfunction. Int. J. Numer. Meth. Fluids, Vol. 48, pp.

集 JST 資料番号:X0587B, Vol. 3, pp. 301–302, 2002.

1023–1040, 2005.

8. 池端昭夫, 正肖鋒. 保存型 cip 法に基づく運動方程式解法および保存型自由界面捕獲スキームを用いた多4. X. Peng F. Xiao, T. Yabe and H.Kobayashi. Conservative and oscillation-less atmospheric transport
6Copyright c 2010 by JSFM

第 24 回数値流体力学シンポジウム
E10-5&"#01234&%"#%$"#$!"#*+,#-./,
*+,#!

%'()

$!

*564789

Fig. 9: Time histories of mean width of interface (viscous ﬂuid)
相流体解析. 日本機械学会計算力学講演会講演論文集

JST 資料番号:L0203A, Vol. 15, pp. 449–450, 2002.
9. 藤岡奨, 牛島省. 非数値拡散フィルタの気液混相流れ数値計算への適用. 水工学論文集 (投稿中).
10. 牛島省, 竹村雅樹, 山田修三, 禰津家久. 非圧縮性流体
解析に基づく粒子 − 流体混合系の計算法 (MICS) の提案. 土木学会論文集, No. 740/II-64, pp. 121–130,2003.
11. 牛島省, 山田修三, 藤岡奨, 禰津家久. 3 次元自由水面
流れによる物体輸送の数値解法 (3D MICS) の提案と適用性の検討. 土木学会論文集, Vol. 810/II-74, pp.
79–89, 2006.

12. Hirsh C. Numerical Computation of Internal andExternal Flows. John Wiley & Sons, 1990.
13. LeVeque R.J. Numerical Methods for Conservation
Laws. Birkhäuser Verlag, 1990.
14. T. Okamoto and M. Kawahara. Two-dimensional
sloshing analysis by the arbitrary lagrangianeulerian ﬁnite element methods. Proceeding ofJSCE, Vol. 441, No. I-18, pp. 39–48, 1992.

7

Copyright c 2010 by JSFM