Transcript 公聴会発表用資料（パワーポイント）

「移流輸送の
高次精度かつ高解像度スキームの開発」
山口大学大学院 理工学研究科
システム工学専攻 社会システム工学大講座
坪郷 浩一
山口大学大学院
坪郷 浩一
＜研究の背景＞
移流計算スキームの開発
数値振動を起こさない
数値拡散が少ない
従来，多くの移流項計算スキームが提案されているがいずれ
も一長一短あり決定版と呼べるものがない．
山口大学大学院
坪郷 浩一
＜目的：移流方程式を数値的に高精度に解く計算手法の開発＞
支配方程式（簡単のため一次元で考える）
対流形式，非保存形式
保存形式
 
u
0
t
x
 u

0
t
x
u=const.の場合に
両式は一致する．
上式を用いて，従来の計算スキームの問題点を指摘する．
Ex1 高次精度スキーム：非保存形式
長所：ガウス分布などの「極値の再現性」に優れている．
短所：「様々な勾配の値を含む分布」の再現性には改良の余地がある．
Ex2 高解像度スキーム：保存形式
長所：「矩形分布での数値振動が抑制」され，再現性に優れている．
短所：ガウス分布などの「極値の再現性」には改良の余地がある．
山口大学大学院
坪郷 浩一
＜研究の目的＞
流れ解析では，基礎方程式中の移流項の計算スキームの離散
化精度が計算解の精度を大きく左右する一つの要因である．本研
究では，数値振動や数値拡散の影響の少ない「高次精度かつ高
解像度スキームの開発」を目指す．
1. Komatsuらが開発した6-point schemeに基づいて高次精度
の6-point schemeを提案するものである．（第2章）
2. 1.で提案した計算スキームを移動境界問題の数値解析に
適用させ有効性の検証を行なう．
（第3章，第4章）
山口大学大学院
坪郷 浩一
＜本研究における高精度，高次精度，高解像度の違い＞
高精度（high-accuracy）：
数値解がより厳密解に近い場合を言う．
高次精度（high-order）：
スキームの打ち切り誤差のオーダーが高いことを意味する
時に用いている．
高解像度（ high-resolution）：
数値振動が発生する箇所を滑らかな単調曲線分布として求
め，それ以外では高精度な解を与えるようなスキーム
いわゆるTVD条件を満足するTVD schemeである．
山口大学大学院
坪郷 浩一
＜各章の関連性＞
移流計算スキームの開発（第1章）
数値振動を起こさない
数値拡散が少ない
高次精度かつ高解像度6-point schemeの開発（第2章）
自由水表面流れ解析
密度関数法（第3章）
Level Set法（第4章）
気液界面の界面を捕捉する移流
計算に第2章で開発した計算スキ
ームを適用させる．
気液界面を精度よく捕捉するために，数値振動や数値拡散が
発生しにくい移流項計算スキームを用いる必要がある．
山口大学大学院
坪郷 浩一
＜第２章 高次精度6-point schemeの開発＞
目的
Komatsuらが開発した1次精度の6-point schemeの導出の概念
に従って，新たな高次精度の6-point schemeの構築を主要な目
的とする．またフラックス制限関数（limiter）が導入できるように保
存形式にも書き換え，さらに朝位らのdiscriminatorの改善を図り，
高次精度の6-point schemeと組み合わせて様々な分布形状に
対し「高次精度かつ高解像度移流計算手法を提案」する．
・特性曲線形式高次精度6-point scheme
・保存形式高次精度6-point scheme
・ 数値振動とクリッピングエラーの除去
・ 多次元問題への適用
山口大学大学院
坪郷 浩一
＜基礎方程式＞
支配方程式（簡単のため一次元で考える）
対流形式，非保存形式
 
u
0
t
x
u=const.の場合に
両式は一致する．
保存形式
山口大学大学院
 u

0
t
x
坪郷 浩一
＜特性曲線法＞
t
ｎ+1
特性曲線
ｎ
x
i-3
 
u
0
t
x
i-2
特性曲線表示

n1
i

n
 
n
i-1
d
0
dt
on
ξ
i+1
i
dx
u
dt
i+2
これは特性曲線上での値
が変わらないことを意味する．
時刻n+1のi点の濃度を求める問題は，特性曲
線によって時刻nの点ξの濃度を内挿によって
求める問題に帰着される．
を近隣の値を用いて内挿式で求める．
山口大学大学院
坪郷 浩一
＜1次精度6-point schemeの最終形＞
in1  P1in3  P2in2  P3in1  P4in  P5in1  P6in2
ここで，
P1  
13 3 3877 2 17117
 
 

720
101280
303840
37 3 1069 2 18821
 
 

144
20256
60768
49
6563 2 31373
P3    3 
 

72
10128
30384
P2 
49 3 4705 2 8717
 
 
 1
72
3376
30384
13 3 5561 2 34435
P5  
 
 

144
6752
60768
P4 
P6  
13 3 3121 2 22603
 
 

720
33760
303840
(=uDt/Dx)はクーラン数である．
山口大学大学院
坪郷 浩一
＜Komatsuらが開発した1次精度の6-point schemeの問題点＞
様々な勾配の値を含む半楕円分布や矩形分
布の再現性には改良の余地がある
ガウス分布などの極値
の再現性に優れている
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
山口大学大学院
坪郷 浩一
＜6 point scheme (Komatsu et al,1985)＞
3次曲線を内挿
x
i-4
i-3
i-2
i-1
Δx
i
i+1
i+2
i+3
時刻nのi-3からi+2までの格子点上の濃度を
i3，i2，i1，i，i1，i2
とする．
境界条件：線形外挿
i 4  2i 3  i 2
i3  2i2 i1
これら8個の濃度の値を使用する．ただし実質には6個の濃度の
値しか使用していないことに注意する．
山口大学大学院
坪郷 浩一
＜3次多項式＞
3次曲線を内挿
x
i-1/2
i-4
i-3
i-2
i-1
i
i+1
i+2
i+3
F ( x)(1)  f  x; i4 , i3 , i2 , i1 
F ( x)(2)  f  x; i3 , i2 , i1, i 
F ( x)(3)  f  x; i2 , i1, i , i1 
F ( x)(4)  f  x; i1, i , i1, i2 
F ( x)(5)  f  x; i , i1, i2 , i3 
4個の濃度値とそのx座標が与えられれば，その濃度値を通過す
る3次多項式F(x)を構築できる．
山口大学大学院
坪郷 浩一
ｉ-1点
ｉ-1点の勾配
ｉ点
ｉ-1/2点の勾配
ｉ点の勾配
x
i-1/2
i-2
i-1
i
i+1
以下のような手順で6-point schemeを構築できる．
   xi 1    xi 1/ 2 

F ( x)  f  x;
,
,   xi 1  ,   xi  
x
x


   xi    xi 1/ 2 

(b )
F ( x)  f  x;
,
,   xi 1  ,   xi  
x
x


(a)
F ( x)(a)  F ( x)(b)
F ( x) 
2
山口大学大学院
6-point schemeの誘導で特
徴的なのはi-1/2の位置の
濃度勾配を使用しているこ
とである．
坪郷 浩一
＜3次精度6-point schemeの導出＞
3次曲線を内挿
x
i-4
i-3
i-2
i-1
i-1/2
i
i+1
i+2
i+3
F ( x)(1)  f  x;  i  4 ,  i 3 ,  i  2 ,  i 1 
F ( x)(2)  f  x;  i 3 ,  i 2 ,  i 1 ,  i 
F ( x)(3)  f  x;  i  2 ,  i 1 ,  i ,  i 1 
F ( x)(4)  f  x;  i 1 ,  i ,  i 1 ,  i  2 
F ( x)(5)  f  x;  i ,  i 1 ,  i  2 ,  i 3 
4個の濃度値とそのx座標が与えられれば，その濃度値を通過す
る3次多項式F(x)を構築できる．
山口大学大学院
坪郷 浩一
＜Φi-1の勾配＞
x
i-1/2
i-3
i-2
i-1
i
i+1
i+2
F ( xi1 )(2)  f   xi1; i3 , i2 , i1, i 
F ( xi1 )(3)  f   xi1; i2 , i1, i , i1 
F ( xi1 )(4)  f   xi1; i1, i , i 1, i 2 
F ( xi1 )   F ( xi1 )(2)   F ( xi1 )(3)  1  2  F ( xi1 )(4)
山口大学大学院
坪郷 浩一
＜Φiの勾配＞
x
i-1/2
i-3
i-2
F ( xi )(2) 
F ( xi )(3) 
i-1
i
i+1
i+2
f   xi ; i3 , i2 , i1, i 
f   xi ; i2 , i1, i , i1 
F ( xi )(4)  f   xi ; i1, i , i1, i2 
F ( xi )  1  2  F ( xi )(2)   F ( xi )(3)   F ( xi )(4)
山口大学大学院
坪郷 浩一
＜Φi-1/2の勾配＞
x
i-1/2
i-3
i-2
i-1
i
i+1
i+2
F ( xi1/ 2 )(3)  f   xi1/ 2 ; i2 , i1, i , i1 
F(xi1/ 2 )  F(xi1/ 2 )(3)
山口大学大学院
坪郷 浩一
in1  P1in3  P2in2  P3in1  P4in  P5in1  P6in2
ここで，
1
1
1
1
P1      2     
6
6
4
4
1
1
2
 3
3
P2   3       2     
6
2
3
 4
4
1
1
4
 1
 1
P3    3       2      
2
6
3
 2
 2
1
4
1
1
 1
P4   3      2       1
2
3
6
2
 2
1
1
5
7
 7
P5    3      2      
6
4
4
 4 12 
1
1
 1
 1
P6       2      
 4 12 
 4 12 
(=uDt/Dx)はクーラン数である．
山口大学大学院
坪郷 浩一
打ち切り誤差解析
Taylor級数展開による打ち切り誤差解析を行えば次式を得る．
なお，打ち切り誤差項は4次までを記す．
0


u
  (1   )(2  1)
t
x
0
0
1  2 Dx2
3  3 Dx3
 1
3
 (1   )      2
  (1   )     3
4  x 2!Dt
4  x 3!Dt
 2
2
 4
 13 17  2  13 13   4 Dx4
3
   2             4
 O(Dx5 )
4
4   x 4!Dt
2
 2

右辺第1項は0次の誤差項である．この項は恒等的に0にな
る必要がある．恒等的に-2 +1=0となれば良いので重み
=1/2が得られる．
山口大学大学院
坪郷 浩一
＜特性曲線形式3次精度6-point scheme＞
in1  P1in3  P2in2  P3in1  P4in  P5in1  P6in2
ここで，
1 2 1
  
24
24
1
1
7
P2   3   2  
6
8
24
1
5
13
P3    3   2  
2
12
12
1
13
5
P4   3   2   1
2
12
12
1
5
11
P5    3   2  
6
8
24
P1  
P6  
1 2 1
  
24
24
(=uDt/Dx)はクーラン数である．
山口大学大学院
坪郷 浩一
＜4次精度6-point schemeの導出＞
3次精度6-point schemeの打ち切り誤差は以下のようになる．
4
4


2   Dx
2
u
  1   
 O(Dx5 )
4
t
x
x 4! Dt
下式を整理すれば，特性曲線形式4次精度6-point schemeである．

n 1
i
 P1
n
i 3
 P2
n
i 2
 P3  P4  P5  P6
n
i 1
n
i
n
i 1
3次精度6-point scheme
山口大学大学院
n
i 2
Dx4 in2  4in1  6in  4in1  in2
  (1   ) 
24
Dx4
2
2
4階微分は中央差分で近似
坪郷 浩一
＜特性曲線形式4次精度6-point scheme＞
in1  P1in3  P2in2  P3in1  P4in  P5in1  P6in2
ここで，
P1  
1 2 1
  
24
24
1 4 1 3 1 2 7
      
24
12
6
24
1
1
1
13
P3    4   3   2  
6
6
4
12
1
5
5
P4   4   2    1
4
6
12
P2 
1
1
11
11
P5    4   3   2  
6
6
24
24
P6 
1 4 1 3 1
    
24
24
24
(=uDt/Dx)はクーラン数である．
山口大学大学院
坪郷 浩一
＜モデル計算＞
計算格子幅 Δx=200m
計算時間間隔 Δt=100sec
流速 0.5m/sec
クーラン数 0.25
初期分布
中心位置1400m，標準偏差264m，ピーク値10
中心位置2400m，標準偏差264m，ピーク値6.5
の二種類のガウス分布の重ね合わせ
幅2000mの矩形分布
中心位置5800mで最大値10の半楕円分布（空間方向の半径2000m）
山口大学大学院
坪郷 浩一
＜モデル計算＞ Δx=200m ， Δｔ= 100sec
12
12
厳密解
計算解
10
8
濃
度
厳密解
計算解
10
8
6
濃
度
4
6
4
2
2
0
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
-2
4000
20000
6000
8000
10000
x軸（ｍ）
14000
16000
18000
20000
x軸（ｍ）
1次精度6-point scheme（9600sec）
3次精度6-point scheme（9600sec）
12
12
厳密解
計算解
10
厳密解
計算解
10
8
8
濃
度
12000
6
濃
度
4
6
4
2
2
0
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
x軸（ｍ）
4次精度6-point scheme（9600sec）
山口大学大学院
20000
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
CIP scheme（9600sec）
坪郷 浩一
＜モデル計算＞ Δx=200m ， Δｔ= 100sec
12
12
厳密解
計算解
10
8
濃
度
8
6
濃
度
4
6
4
2
2
0
0
-2
95200
厳密解
計算解
10
97200
99200
101200
103200
105200
107200
109200
111200
-2
95200
97200
99200
101200
x軸（ｍ）
107200
109200
111200
3次精度6-point scheme（192000sec）
12
12
厳密解
計算解
10
厳密解
計算解
10
8
8
6
濃
度
4
6
4
2
2
0
0
-2
95200
105200
x軸（ｍ）
1次精度6-point scheme（192000sec）
濃
度
103200
97200
99200
101200
103200
105200
107200
109200
x軸（ｍ）
4次精度6-point scheme（192000sec）
山口大学大学院
111200
-2
95200
97200
99200
101200
103200
105200
107200
109200
111200
x軸（ｍ）
CIP scheme（192000sec）
坪郷 浩一
＜モデル計算＞
計算格子幅 Δx=100m クーラン数一定のまま空間解像度を倍にする．
計算時間間隔 Δt=100sec
流速 0.25m/sec
クーラン数 0.25
初期分布
中心位置1400m，標準偏差264m，ピーク値10
中心位置2400m，標準偏差264m，ピーク値6.5
の二種類のガウス分布の重ね合わせ
幅2000mの矩形分布
中心位置5800mで最大値10の半楕円分布（空間方向の半径2000m）
山口大学大学院
坪郷 浩一
＜モデル計算＞ Δx=100m ， Δｔ= 100sec
12
12
厳密解
計算解
10
8
8
濃
度
厳密解
計算解
10
6
濃
度
4
6
4
2
2
0
0
-2
95200
97200
99200
101200
103200
105200
107200
109200
-2
95200
111200
97200
99200
101200
1次精度6-point scheme（192000sec）
105200
107200
109200
111200
3次精度6-point scheme（192000sec）
12
12
厳密解
計算解
10
8
濃
度
6
4
厳密解
計算解
10
8
濃
度
103200
x軸（ｍ）
x軸（ｍ）
6
4
2
2
0
0
-2
95200
97200
99200
101200
103200
105200
107200
109200
x軸（ｍ）
4次精度6-point scheme（192000sec）
山口大学大学院
111200
-2
95200
97200
99200
101200
103200
105200
107200
109200
111200
x軸（ｍ）
CIP scheme（192000sec）
坪郷 浩一
＜モデル計算＞ Δx=100m ， Δｔ= 100sec
12
12
厳密解
計算解
10
8
8
濃
度
厳密解
計算解
10
6
濃
度
4
6
4
2
2
0
0
-2
95200
97200
99200
101200
103200
105200
107200
109200
-2
95200
111200
97200
99200
101200
1次精度6-point scheme（192000sec）
105200
107200
109200
111200
3次精度6-point scheme（192000sec）
12
12
厳密解
計算解
10
8
濃
度
6
4
厳密解
計算解
10
8
濃
度
103200
x軸（ｍ）
x軸（ｍ）
6
4
2
2
0
0
-2
95200
97200
99200
101200
103200
105200
107200
109200
x軸（ｍ）
4次精度6-point scheme（192000sec）
山口大学大学院
111200
-2
95200
97200
99200
101200
103200
105200
107200
109200
111200
x軸（ｍ）
CIP scheme（192000sec）
坪郷 浩一
＜モデル計算＞ Δx=100m ， Δｔ= 100sec
ガウス分布の最大値8.585
12
ガウス分布の最大値7.935
厳密解
計算解
10
12
8
8
濃
度
厳密解
計算解
10
6
濃
度
4
6
4
2
2
0
0
-2
95200
97200
99200
101200
103200
105200
107200
109200
x軸（ｍ）
4次精度6-point scheme（192000sec）
111200
-2
95200
97200
99200
101200
103200
105200
107200
x軸（ｍ）
CIP scheme（192000sec）
109200
111200
＜保存形式高次精度6-point scheme ＞
矩形分布で数値振動が発生している．
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
特性曲線形式高次精度6-point schemeにフラックス制限関数を導入する．
なお，これ以後の特性曲線形式高次精度6-point schemeは，TVD化するため
に保存形式に直すため保存形式高次精度6-point schemeと呼ぶことにする．
CIP schemeは保存形式に直すことができず，TVD化できない．
山口大学大学院
坪郷 浩一
＜保存形式高次精度6-point schemeへの拡張 ＞
保存形式の移流方程式を有限体積法で離散化
 u

0
t
x
離散化
整理すれば


n
i
n
i 1/ 2
    (
n
i
ui-1/2
u ：流速一定
 u

Dt
n 1
i
n 1
i
コント ロールボリューム
n
i 1/ 2
 u
Dx
n
i 1/ 2

n
i 1/ 2
0
i-1
ui+1/2
i
i - 1/2
左側界面
)
ｘ
i+1
i + 1/2
右側界面
n
i
1/ 2 はコントロールボリュームの界面濃度だから5点の濃度で以下のように近似できる．
in1/ 2  C1in2  C2in1  C3in  C4in1  C5in2
in1/ 2  C1in3  C2in2  C3in1  C4in  C5in1
まとめると
in1  C1in3  (C1  C2 )in2  (C2  C3 )in1  (C3  C4 ) 1in  (C4  C5 )in1 C5in2
(=uDt/Dx)はクーラン数である．
山口大学大学院
坪郷 浩一
＜本論文で提案する保存形式高次精度6-point scheme＞
in1  C1in3  (C1  C2 )in2  (C2  C3 )in1  (C3  C4 ) 1in  (C4  C5 )in1 C5in2
1次精度6-point scheme
13 2
3877
17117
C1  
 

720
101280
303840
43 2
367
19247
C2 
 

180
25320
75960
53 2 53
187
C3  
  
120
80
240
43 2 6171
37393
C4 
 

180
8440
75960
13 2 3121
22603
C5  
 

720
33760
303840
3次精度6-point scheme
1
1
C1    
24
24
1 2 1
1
C2     
6
12
4
1 2 1
5
C3      
3
2
6
1 2 7
5
C4     
6
12
12
1
1
C5   
24
24
4次精度6-point scheme
1
1

24
24
1
1
1
1
 3  2   
24
12
8
4
1 3 1 2 3
5
     
8
12
8
6
1 3 1 2 11
5
     
8
12
24
12
1
1
1
  3   
24
12
24
C1  
C2
C3
C4
C5
(=uDt/Dx)はクーラン数である．
山口大学大学院
坪郷 浩一
＜モデル計算＞
計算格子幅 Δx=200m
計算時間間隔 Δt=100sec
流速 0.5m/sec
クーラン数 0.25
初期分布
中心位置1400m，標準偏差264m，ピーク値10
中心位置2400m，標準偏差264m，ピーク値6.5
の二種類のガウス分布の重ね合わせ
幅2000mの矩形分布
中心位置5800mで最大値10の半楕円分布（空間方向の半径2000m）
山口大学大学院
坪郷 浩一
＜universal limiterの導入＞ Δx=200m ， Δｔ= 100sec
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
特性曲線形式1次精度6-point scheme
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
保存形式1次精度6-point scheme
universal limiter適用
山口大学大学院
坪郷 浩一
＜universal limiterの導入＞ Δx=200m ， Δｔ= 100sec
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
特性曲線形式1次精度6-point scheme
12
厳密解
計算解
10
8
濃
度
6
4
数値振動は除去されているが，ガウスの
極値の再現性の低下が問題
なお，極値の再現性の低下をクリッピン
グエラーと呼ぶ．
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
x軸（ｍ）
保存形式1次精度6-point scheme
universal limiter適用
山口大学大学院
20000
数値振動だけを除去でき，物理的に意味の
ある極値近傍ではlimiterを機能させない朝
位ら(１９９９）のdiscriminatorを導入する．
坪郷 浩一
＜朝位らのdiscriminatorの導入＞
流速の方向
D1×D3 ＜0
| D2 |＜| D1 |
No
limiter
on
極値
実際の分布
Yes
D3
limiter off at i+1/2, i-1/2
  i2
   i 1
  i
D1  i 1
D2  i
D3  i 1
Dx
Dx
Dx
D2
D1
i 1/ 2
朝位らのdiscriminatorのアルゴリズム
i2
i 1
i  1/ 2
limiter off
i
i 1
朝位らのdiscriminatorの概念図
山口大学大学院
坪郷 浩一
＜朝位らのdiscriminatorの導入＞ Δx=200m ， Δｔ= 100sec
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
保存形式1次精度6-point scheme
universal limiter適用
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
保存形式1次精度6-point scheme
universal limiterとdiscriminatorの適用
山口大学大学院
坪郷 浩一
＜朝位らのdiscriminatorの導入＞ Δx=200m ， Δｔ= 100sec
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
ガウスの極値の再現性は向上したが矩形分
布の裾の近傍で数値振動が発生している．
保存形式1次精度6-point scheme
universal limiter適用
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
保存形式1次精度6-point scheme
universal limiterとdiscriminatorの適用
山口大学大学院
坪郷 浩一
＜朝位らのdiscriminatorの誤作動＞
D1×D3 ＜0
| D2 |＜| D1 |
No
limiter
on
Yes
D3  0
limiter off at i+1/2, i-1/2
D1 
 i 1   i  2
   i 1
  i
D2  i
D3  i 1
Dx
Dx
Dx
急勾配
D≒　1 0
 D1  0
i-2
D2  0
i-1
i
i+1
朝位らのdiscriminatorが「極値がある」と誤作動したために数値振動が発生した．
改善
「勾配D1とD3が大きく異なる場合」においてlimiterを作用させない．
山口大学大学院
坪郷 浩一
＜修正された朝位らのdiscriminator＞
D1×D3＜0
| D2|＜| D1|
No
limiter
on
Yes
Max  D1 ,  
Max  D3 ,  

Yes

Yes
Maxは2個の引数の中で大きい方を選択す
る演算子である．は分母が0にならないよ
うにするために導入される非常に小さな値
である．は任意の閾値である．
No
Max  D3 ,  
Max  D1 ,  
No
limiter off at i+1/2, i-1/2
修正された朝位らの
discriminatorのアルゴリズム
山口大学大学院
坪郷 浩一
＜修正された朝位らのdiscriminator＞
移流問題を用いて提案した修正型discriminatorの機能の評価を
保存形式1次精度6-point schemeを対象に閾値の検討を行った．
ガ
ウ
ス
分
布
の
ピ
ー
ク
値
10
0.00
9
計
算
領
域
内
の
最
小
値
Limiter を作用させない
場合のピーク値
8
Limiter を作用させた
場合のピーク値
7
6
0
5
10
閾値ε
15
20
閾値とガウス分布のピーク値の関係
-0.05
-0.10
-0.15
-0.20
limiter を作 用さ せ
-0.25
ない場合の最小値
-0.30
-0.35
0
5
10
15
20
閾値ε
閾値と最小値の関係
数値実験によりは，3≦ 1が推奨される．本論文では， =3を用いる．
山口大学大学院
坪郷 浩一
＜修正された朝位らのdiscriminator＞
12
12
厳密解
計算解
10
厳密解
計算解
10
8
濃
度
Δx=200m ， Δｔ= 100sec
8
6
6
濃
度
4
4
2
2
0
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
-2
4000
6000
8000
10000
x軸（ｍ）
14000
16000
18000
20000
x軸（ｍ）
特性曲線形式1次精度6-point scheme
保存形式1次精度6-point scheme
universal limiter適用
12
12
厳密解
計算解
10
厳密解
計算解
10
8
8
濃
度
12000
6
濃
度
4
6
4
2
2
0
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
x軸（ｍ）
保存形式1次精度6-point scheme
universal limiterとdiscriminatorの適用
山口大学大学院
20000
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
保存形式1次精度6-point scheme
universal limiterと修正型discriminatorの適用
坪郷 浩一
＜保存形式高次精度6-point scheme ＞ universal limiterと修正型discriminatorの適用
12
12
厳密解
計算解
10
8
濃
度
8
6
濃
度
4
6
4
2
2
0
0
-2
4000
厳密解
計算解
10
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
保存形式1次精度6-point scheme
保存形式3次精度6-point scheme
12
厳密解
計算解
10
8
濃
度
6
4
2
0
-2
4000
6000
8000
10000
12000
14000
16000
18000
20000
x軸（ｍ）
保存形式4次精度6-point scheme
山口大学大学院
坪郷 浩一
＜保存性の検証＞

(t ) 
i i max
VErr
i 1

Ci (t )  Ci (0)
i i max
i 1
100
Ci (0)
ここでCi (t)は計算終了時刻t におけるのi点の濃度，Ci (0)は初期
のi 点の濃度，imax は計算点の総数である．
質量の相対誤差は，10-6%程度である，limiterやdiscriminatorを導入して
も保存性は良好である．
山口大学大学院
坪郷 浩一
＜多次元問題＞
～剛体回転する流れ場におけるモデル計算～
2次元平面で角速度ω=2π/12000（sec-1）
計算格子幅 Δx= Δｙ＝１00m
計算時間間隔 Δt=50sec
計算領域 ４０００×４０００
初期分布
中心位置（1400m，1400m） 標準偏差200m，ピーク値10
のガウス型濃度分布
計算スキーム
特性曲線形式4次精度6-point scheme
保存形式4次精度6-point scheme
QUICKEST scheme
limiterやdiscriminatorを適用した．
山口大学大学院
坪郷 浩一
～剛体回転する流れ場における移流計算の結果～ (x=y上
で)
ピーク値の再現性は特性曲線6-point schemeが最も良い．
12
10
厳密解
保存形式
特性曲線型
QUICKEST
8
濃
度
6
4
2
0
0
1000
2000
3000
4000
-2
距離（m）
質量の相対誤差は保存形式6-point schemeが最も良い．
質量の相対誤差
保存形式4次精度6-point scheme：3.95×10-3%
特性曲線形式4次精度6-point scheme ： 7.56×10-2%
QUICKEST scheme：3.43×10-2%
山口大学大学院
坪郷 浩一
＜高次精度6-point schemeの非線形移流項＞
～2次元キャビティ問題～
Ghia et al.の結果と比較できるように正方キャビティを計算対
象とした．MAC法の一種であるFractional Step法で数値計算
を行った．
計算条件
等間隔で50×50のセルに分割しstaggered格子を用いる．
Re=1000
計算スキーム
保存形式3次精度6-point scheme
保存形式4次精度6-point scheme
5次精度風上差分法
limiterを適用した．
山口大学大学院
坪郷 浩一
～2次元キャビティ内の流動計算結果～
キャビティ中心を含む鉛直方向（y方向)および水平方向(x方向)の流速分布
1.0
1.0
0.5
x
ｙ
Ghia
４次精度6point法
５次精度風上差分法
0.5
Ghia
4次精度6point法
5次精度風上差分法
0.0
-0.4
-0.2
0.0
0.2
0.4
u
x方向流速uのy方向分布
山口大学大学院
0.6
0.0
-0.5
0.0
0.5
v
y方向流速vのx方向分布
坪郷 浩一
高次精度かつ高解像度6-point schemeの開発（第2章）
本章では，高次精度かつ高解像度6-point schemeの開発を行い，以下のような知見
が得られた．
1. 極値の再現能力は，CIP scheme より優れている．
2.フラックス制限関数と修正されたdiscriminatorを用いれば，高精度かつ高解像度数値
解を得ることが可能となった．一方，CIP schemeは保存形式に直すことができず，TVD
化できない．
3. 高次精度6-point scheme はtime splitting法を用いることで 多次元問題にも容易に適
用でき，精度の良い解が得られることが示された．
2.の操作が人為的な側面を持っている．今後の課題として開発した計算ス
キームのみで単調性や凹凸性の保存ができるような計算スキームの開発を
検討する必要がある．
山口大学大学院
坪郷 浩一
高次精度かつ高解像度6-point schemeの開発（第2章）
＜適用事例（第3章，第4章）＞
自由水表面流れ解析
気液界面を捕捉する関数の移流項に，提案した計算スキーム
を使用する．
界面捕捉法の有力な手法の
密度関数法 （第3章）
Level Seｔ法（第4章）
を用いて提案した計算スキームの有効性を移動境界流れの検証
モデルの1つであるダムブレイクモデルで検証する．
山口大学大学院
坪郷 浩一
＜第3章 密度関数法を用いた自由水表面流れ解析＞
目的
本章では，密度関数の移流方程式に提案した高次精度差分
スキーム（第2章）を適用させその有効性の検証を行う．また，
計算進行に伴う界面のぼやけを回避し体積を補正する手法の
開発について述べるものである．
・ 数値計算法
液相体積の計算
体積補正法
質量補正法
・ 提案した計算スキーム（第2章）の検証
・ まとめ
山口大学大学院
坪郷 浩一
＜自由水表面が存在する流れの数値解析手法＞
界面捕捉法
複雑で大規模な自由水表面問題に対しては界面捕捉法が有利である．
代表例
・VOF法(1981，C.W.Hirtら）：液相のみ
最近の研究：気液二相流を直接解く手法が発展している
・密度関数法（1996，Kanaiら）
・Level Set法（1994，Sussmannら）
山口大学大学院
坪郷 浩一
界面捕捉法の問題点：体積の保存性に問題がある．
既往の体積補正の手法として，
Level Set法では姫野らの手法（1999年）
VOF法では桜庭らの手法（2003年）
が挙げられるが，密度関数法に関する体積補正の手法は，まだ
開発されていない．
山口大学大学院
坪郷 浩一
＜ 目 的 ＞
本研究では，自由水表面を持つ流れ解析の有力な手法の一つ
である密度関数法を用いて以下のような解析・検討を行った．
＜移流項計算スキームの検証＞
本章では，密度関数の移流方程式に提案した高次精度差分ス
キーム（第2章）を適用させその有効性の検証を行う．
＜体積補正の手法の提案＞
① この手法のための「オリジナルの体積補正法」を開発する．
② ①で必要な液相体積を求めるアルゴリズムを開発する．
山口大学大学院
坪郷 浩一
＜基礎方程式＞
連続の式
 u  0
運動方程式
u
1
 u  u   F p  2u
t

密度の保存式

  u   0
t
ここでuは流速ベクトル，pは圧力， μは粘性係数，Fは体積力である．
山口大学大学院
坪郷 浩一
＜密度の保存式＞

  u   0
t
密度の保存式を直接解く代わりに密度関数（０＜φ＜１）の保存則を解く

  u   0
t
密度関数Φ（0<φ<1）の判別
気相Φ=0，界面Φ=0.5，液相Φ=1
密度および粘性係数
  Liq  (1 )Gas
  Liq  (1 )Gas
ここに ρLip は液相の密度 ，ρGasは気相の密度 ，
μLip は液相の粘性係数 ，μGas は気相の粘性係数である．
山口大学大学院
坪郷 浩一
＜体積補正法の概念＞
界面セルの液相占有率が真の
それよりも小さいために液相体
積が過小評価する
界面の総体積
真の界面
界
面セルの総体積 真
の界面
界面セルの液相占有率を
高くする必要がある．
V
Err
液相体積誤差量V
Eｒｒ 計計算界面
算界面
液相占有率の補正量を算出し，
それを各界面セルに加えること
で体積誤差を減少させる．
山口大学大学院
坪郷 浩一
＜界面セルの把握状態の一例＞
気液界面セル
気相セル
液相セル
Φ=0.5の等高線
＜液相体積（面積）を求める＞
VLiq (t)   D()dxdy
山口大学大学院
1
液相セル


D()    x, y, t  界面セル

0
気相セル

坪郷 浩一
＜体積補正の手法＞
現在時刻ｔの液相体積誤差量VErr（ｔ）を求める
VErr  VLiq (t ) -Vini
VLiq(t)：任意時刻tの液相体積
Vini：初期液相体積
体積補正量LErrを求める
LErr 
VErr (t )
A(t )
A(t)：任意事項tの計算領域内の界面セルの総体積
Φn+1を求める時点で得られた密度関数に次のような操作を行う．

n1
山口大学大学院

n1
Lerr  0 界面･気相セル
 Lerr 
Lerr  0 界面･液相セル
坪郷 浩一
＜質量補正の導入＞
体積補正法の導入により「計算領域内の全質量保存を破綻」
させてしまう可能性がある．
原因
質量誤差は界面・界面近傍セルのΦ値が適切でないために生じる
ものと考えられる．
改善点
「体積補正の後に全質量の保存」を保証する目的のルーチンを作成
山口大学大学院
坪郷 浩一
＜質量補正の手法＞
現在時刻ｔの全質量誤差量MErr（ｔ）を求める
MErr (t)  M (t)  Mini
M(t)：任意時刻tの全質量
Mini：初期全質量
質量補正量MErrを求める．
M Err (t )
M Err 
AM (t )
MErr(t)：任意時刻tの質量誤差
AM(t)：界面近傍の0＜Φ＜1となる
領域の任意時刻ｔの体積
計算領域内のセルから質量誤差を一様に差し引く．

山口大学大学院
n1

n1
 M Err
坪郷 浩一
＜計算のアルゴリズム＞
START
仮の流速u*，ｖ*を求める
圧力のポアソン方程式を解き
n+1
ｐ を求める
u
n+1
ｎ+1
，ｖ
を求める
Φｎ+1を求める
体積補正・質量補正を行う
体積補正法
ρ
NO
n+1
，μ
ｎ+1
を求める
計算時間の終了時刻
に達したか？（t>nΔt）
YES
END
山口大学大学院
坪郷 浩一
＜移流項計算スキームの検証＞
空気
密度：1.25kg/㎥
粘性係数：1.50×10-5Pa・ｓ
水
密度：1000kg/㎥
粘性係数：1.00×10-3Pa・ｓ
Ｌ＝0.15ｍ，0.6ｍ×0.6m
計算領域40×40
計算格子間隔Δx=Δy = 0.015m
時間刻みは Δt = 0.0001sec
表面張力なし
山口大学大学院
坪郷 浩一
＜移流項計算スキームの検証＞
運動方程式
TVD
輸送方程式
TVD
CIP
安定化CIP
3C6P
TVD
CIP
安定化CIP
3C6P
山口大学大学院
体積補正法なしで最も数値拡散が少
ない組み合わせの検討
：3次精度TVD-MUSCL scheme
：CIP scheme
：安定化CIP scheme
：3次精度保存形式6-point scheme
坪郷 浩一
＜全質量の相対誤差＞
全質量の相対誤差(％)
10
5
TVD-TVD
TVD-3C6P
TVD-安定化CIP
0
0
山口大学大学院
2
4
6
8
10
時間t(s)
12
14
16
18
20
坪郷 浩一
液相体積の相対誤差(％)
＜液相体積の相対誤差＞
100
TVD-TVD
TVD-3C6P
TVD-安定化CIP
80
60
40
20
0
0
2
4
6
8
10
12
14
16
18
20
時間t(s)
液相体積・全質量の相対誤差からTVD-3C6Pの
組み合わせが最も数値拡散が少ない．
山口大学大学院
坪郷 浩一
＜液相挙動 静止画＞ TVD-3C6P
0.6
0.6
0.6
0.6
0.3
0.3
0.3
  0.001
0.3
Φ = 0.5
0
0
0.3
0.6
0
0
t=0.10s
0.3
0.6
0
0
t=0.30s
0.3
0.6
0
0
t=3.0s
0.3
0.6
t=5.0s
(a)体積補正法なし
0. 6
0.6
0. 3
0.3
0
0
0. 3
t=0.10s
0. 6
0
0
Φ0 .=6 0.5
0. 6
0.3
0.3
0.6
t=0.30s
0
0. 3
0
0.3
t=3.0s
0.6
0
0
0. 3
0. 6
t=5.0s
(b)体積補正法あり
山口大学大学院
坪郷 浩一
＜液相挙動 動画＞ TVD-3C6P
体積補正法なし
山口大学大学院
体積補正法あり
坪郷 浩一
4
4
3
3
Z/L
Z/L
＜水柱先端位置＞
2
Exp.(Martin & Moyce, 1.25in)
Exp.(Martin & Moyce, 2.25in)
Koshizuka et al
1
2
Exp.(Martin & Moyce, 1.25in)
Exp.(Martin & Moyce, 2.25in)
Koshizuka et al
1
体積補正法なし
体積補正法あり
体積補正法なし
体積補正法あり
0
0
0
2
ｔ(2g/L)1/2
TVD-TVD
山口大学大学院
4
0
2
ｔ(2g/L)1/2
4
TVD-3C6P
坪郷 浩一
＜液相体積と質量の相対誤差＞ TVD-3C6P
全質量の相対誤差(％)
10
10
10
10
10
10
10
10
10
10
-5
全質量の相対誤差
-6
-6
-6
-6
-6
-6
-6
-6
-6
0
100
0
2
4
6
8
10
12
14
16
18
20
液相体積の相対誤差(％)
時間t(s)
10
10
10
10
10
10
10
液相体積の相対誤差
-4
-4
-4
-5
-5
-5
-5
0
100
0
2
山口大学大学院
4
6
8
10
時間t(s)
12
14
16
18
20
坪郷 浩一
＜第3章のまとめ＞
第2章で提案した計算スキームの有効性の検証結果を以下に示す．
体積補正法を施さない場合
界面を補足する移流計算に提案した計算スキーム（3次精度保
存形式6-point scheme）を適用したところ液相体積の保存能力
が他の計算スキームと比べて高い，それでも十分な保存性が得
られない．
体積補正法を施した場合
本研究で提案する体積補正法の適用により体積保存と質量保
存の精度が向上することを確認した．しかしながら，流動が抑制
される傾向も現れる．この問題は密度関数の輸送方程式に計算
精度の良い提案した計算スキーム（第2章）を用いることである程
度回避することができる．
山口大学大学院
坪郷 浩一
＜まとめ＞
1. 極値の再現能力は，CIP scheme より優れている．
2. フラックス制限関数と修正されたdiscriminatorを用いれば，高精度かつ高解像度数値
解を得ることが可能となった．
3. 高次精度6-point scheme はtime splitting法を用いることで 多次元問題にも容易に適
用でき，精度の良い解が得られることが示された．
4. 自由水表面流れ解析の一つである密度関数法の輸送方程式の移流項に提案した計
算スキームの適用事例を以下に示す．なお，モデル計算にはダム破壊現象を用いた．
（1）体積補正法を適用しない場合
従来の計算スキーム比べて界面近傍のぼやけや液相体積の減少をある程度抑制し
ている．
（2）体積補正法を適用した場合
従来の計算スキームでは，流体運動が抑制される傾向があった．しかし，提案した計
算スキームを用いることで流体運動の抑制をある程度回避することができた．
山口大学大学院
坪郷 浩一