第 27 回数値流体力学シンポジウム
E04-3

非構造直交格子における空間 5 次精度補間法の構築
Construction of 5th-Order Accurate Spatial Interpolation Method
on Unstructured Cartesian Grids○ 玉置義治, 東大院, 東京都文京区本郷 7-3-1, 4644981114@mail.ecc.u-tokyo.ac.jp
今村太郎, 東大工, 東京都文京区本郷 7-3-1, timat@mail.ecc.u-tokyo.ac.jp
Yoshiharu Tamaki, University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, 113-8656, Japan
Taro Imamura, University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, 113-8656, Japan
In this study, a new high-order interpolation method for unstructured Cartesian grids is developed. This scheme istypically for finite-volume method (FVM). The gradients conventionally calculated by Green-Gauss method arereformed with the value of only the adjacent cells to the target, and then 5th-order accurate interpolated value on thecell surfaces is achieved in even interval grids. This scheme is also adaptive to non-even interval cells. One- andtwo-dimensional verification problems are conducted and it is proved that this scheme is less dissipative than theconventional methods. Also, no significant unphysical reflection or dissipation is observed at the hanging-nodes.
Additionally, it is notable that the computation time is only 1.1 times larger compared to the conventional 3rd-order
upwind-biased scheme.
１．序論
流体の数値計算を行うためには，従来は構造格子を用いるのが一般的であった．
構造格子では各方向のセルの順番(2 次元では i, j)が明示されており，2 つ以上離れたセルの情報を用いて高次精度の計算スキームを構築するのが容易である．しかし，複雑な形状
周りに格子を作成する際には，格子の歪みが問題となり，また格子生成にかかる時間も増大する．
上記の問題点を解決するため，近年四分木法を用いた直交格子の自動生成が注目されており，当研究室においてもソルバーの開発が行われてきた(1)(2)．このソルバーによって自動生成された格子(Fig.1 参照)は，格子の角がすべて直角であるため，ヤコビアンや
メトリックスの評価が不必要であるという利点がある．また解像度の必要性に応じて，任意の箇所に Refinement Box を作成することもできる．一方でハンギングノードがあるため，一般的には非構造格子として扱われ有限体積法による計算が行われる．
非構造格子ではステンシルを伸ばした先に 2 つ隣のセルが存在するとは限らないため，セルに格納されている基本変数(密度，流速，圧力)の勾配を各セルが持つことにより，セル表面での基本変数の値を 2 次精度以上で内挿し，近似 Riemann ソルバーを用いて計算を行う．ただし，勾配計算時には隣接するセルのみしか参照できないため，一般的に勾配値は空間 2 次精度(内挿値は空間 3 次
精度)が限界とされる(3), (4)．しかし，音響等流体の非定常現象を解析する際には，空間 2 次精度や 3 次精度では散逸が大きく，より高次精度のスキームを用いることで大幅な格子数の削減が可能となる．また，Discontinuous Galerkin (DG)法(5)や Spectral Volume 法(6)による非構造格子の高次精度化も試みられているが，演算量の増大が問題となる．
本研究では，直交格子の利点を生かし，勾配値を高次精度で求めることによりセル表面での内挿値を高次精度化する方法を構築する．物体近傍や遠方に存在するハンギングノードにおいては，精度は幾分低下したとしても同じスキームで破綻なく計算が行えるようにする．さらに，この手法を用いて 1 次元，2 次元の問題
において，精度の検証および非物理的な現象が生じないかどうかの確認を行う．また，計算時間を測定し，従来の計算スキームを用いた場合と比較を行う．

Fig.1 Typical example of unstructured Cartesian grid

２． 計算手法
２．１．有限体積法における補間スキームについて
有限体積法では積分型の支配方程式(式(1))が用いられる．
∂∫ 𝑄𝑑𝑉 + ∫ 𝐅 ∙ 𝐧 𝑑𝑆 = 0
∂𝑡 𝑉
𝑆

(1)

第 1 項は保存量𝑄のセル内空間𝑉内での積分を，第 2 項は流束𝐅のセル表面𝑆上での面積分を表している．各セルは変数のセル平均値をデータとして持っているので，単純に体積を掛けることで第 1 項の体積分が行われる．第 2 項は，セル平均値からセル表面に内挿を行い，セル表面両側の値から Riemann 流束を評価して，面積を掛けることで積分される．ここでの空間精度は，(1)内挿の精度，(2)面積分の精度，の 2 つから決定される．前者はセルが持つ勾配量の評価の精度，内挿式による．また，後者については通常，面中心の値に面積を掛けるが，これは空間 2 次精度に留まる．
ここではまず， 1 次元問題において前者の誤差のみを取り扱う．
変数のセル平均値はセル中心での変数値とは異なるため，セル表面へ内挿を行う際には注意が必要となる．ここでは，Fig.2 のよ
うな1次元等間隔の格子においてセル内でのある変数𝑢(𝑥)の分布を考える．以下，j 番目のセル中心での変数値を𝑢𝑗 ，𝑛階の勾配値を∇𝑛 𝑢𝑗 とし，セル平均の値には上線を付けて𝑢𝑗 ，∇𝑛 𝑢𝑗 のように表す．
1

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E04-3
𝑢𝑗+1 − 𝑢𝑗−1
(11)
∇𝑢𝑗 =
+ 𝑜(Δ𝑥 2 )Δ𝑥𝑢𝑗+1 − 2𝑢𝑗 + 𝑢𝑗−1
(12)
∇2 𝑢𝑗 =+ 𝑜(Δ𝑥 2 )Δ𝑥 2としている．ただしこれは等間隔格子における定式化であり，不等間隔部において完全な空間 3 次精度とはなっていないことに注意が必要である．

Fig.2 Even interval one-dimensional gridまず，𝑢(𝑥)を点 j 周りに Taylor 展開すると，
12𝑢(𝑥) = 𝑢𝑗 + ∇𝑢𝑗 (𝑥 − 𝑥𝑗 ) + ∇2 𝑢𝑗 (𝑥 − 𝑥𝑗 )
2
1
3
(2)
+ ∇3 𝑢𝑗 (𝑥 − 𝑥𝑗 )61
4
+ ∇4 𝑢𝑗 (𝑥 − 𝑥𝑗 ) + 𝑜(Δ𝑥 5 )24となる．式(2)を j-1/2～j+1/2 の範囲で積分すると，変数𝑢のセル平均値は式(3)のように求められる．

２．２．空間 5 次精度補間スキームの構築
ここでは等間隔格子を用いることを前提とし，空間 5 次精度の補間スキームを導出する．なお，不等間隔部における対処については後で述べる．
式(2)を積分すると，
j+1 およびj+2 番目のセルでの値は，式(13)，(14)のように計算される．
𝑢𝑗+1 = ∫

3
𝑗+
2

1
𝑗+
2

2

1
1
= 𝑢𝑗 + ∇2 𝑢𝑗 Δ𝑥 2 +∇4 𝑢𝑗 Δ𝑥 4 + 𝑜(Δ𝑥 6 )24
1920
同様にして，勾配のセル平均値は以下のように求まる．
1∇𝑢𝑗 = ∇𝑢𝑗 + ∇3 𝑢𝑗 Δ𝑥 2 + 𝑜(Δ𝑥 4 )24
1 4∇2 𝑢𝑗 = ∇2 𝑢𝑗 +∇ 𝑢𝑗 Δ𝑥 2 + 𝑜(Δ𝑥 4 )24

13 2∇ 𝑢𝑗 Δ𝑥 2
24

5 3
121 4
∇ 𝑢𝑗 Δ𝑥 3 +∇ 𝑢𝑗 Δ𝑥 4
24
1920+ 𝑜(Δ𝑥 5 )

1

1 𝑗+2
𝑢𝑗 =
∫ 𝑢(𝑥)𝑑𝑥
Δ𝑥 𝑗−1𝑢(𝑥)𝑑𝑥 = 𝑢𝑗 + ∇𝑢𝑗 Δ𝑥 +

(13)

+
(3)
𝑢𝑗+2 = ∫

5
𝑗+
2

3
𝑗+
2

(5)
(6)

∇4 𝑢𝑗 = ∇4 𝑢𝑗 + 𝑜(Δ𝑥 2 )

(7)

49 2
∇ 𝑢𝑗 Δ𝑥 2
24

(14)
17 3
1441 43
4+ ∇ 𝑢𝑗 Δ𝑥 +∇ 𝑢𝑗 Δ𝑥12
1920
5+ 𝑜(Δ𝑥 )さらに同様にして j-1，j-2 番目のセル平均値も算出して用いると，最終的に 1，2，4 階の微分値は式(15)～(17)のように表される(1 次元ではセル中点の値を用いたときと等価)．
4 𝑢𝑗+1 − 𝑢𝑗−1
1 𝑢𝑗+2 − 𝑢𝑗−2
(15)
∇𝑢𝑗 = ()− () + 𝑜(Δ𝑥 4 )
3
2Δ𝑥
3
4Δ𝑥

(4)

∇3 𝑢𝑗 = ∇3 𝑢𝑗 + 𝑜(Δ𝑥 2 )𝑢(𝑥)𝑑𝑥 = 𝑢𝑗 + 2∇𝑢𝑗 Δ𝑥 +式(4)～(7)に示された j 番目のセルにおけるセル平均値を用いてセル内の変数𝑢(𝑥)の分布を求めると，最終的に式(2)は式(8)のように整理できる．
1Δ𝑥 2
2𝑢(𝑥) = 𝑢𝑗 + ∇𝑢𝑗 (𝑥 − 𝑥𝑗 ) + ∇2 𝑢𝑗 {(𝑥 − 𝑥𝑗 ) −
}
2
12
1
Δ𝑥 2
3
+ ∇3 𝑢𝑗 {(𝑥 − 𝑥𝑗 ) −(𝑥 − 𝑥𝑗 )}
6
4
(8)
2
1 4
Δ𝑥
7Δ𝑥 4
4
2
+ ∇ 𝑢𝑗 {(𝑥 − 𝑥𝑗 ) −(𝑥 − 𝑥𝑗 ) +
}
24
2
240
+𝑜(Δ𝑥 5 )ここで 𝑥 = 𝑥𝑗 + Δ𝑥⁄2を代入すると，セル表面での値は式(9)のように表される．
1
1
1 4
𝑢 1 = 𝑢𝑗 + ∇𝑢𝑗 Δ𝑥 + ∇2 𝑢𝑗 Δ𝑥 2 −∇ 𝑢𝑗 Δ𝑥 4
𝑗+
2
12
720
(9)
2
+ 𝑜(Δ𝑥 5 )4 𝑢𝑗+1 − 2𝑢𝑗 + 𝑢𝑗−1
1 𝑢𝑗+2 − 2𝑢𝑗 + 𝑢𝑗−2
∇2 𝑢𝑗 = (
)− (
)
3
Δ𝑥 2
3
4Δ𝑥 2
+𝑜(Δ𝑥 4 )(16)

𝑢𝑗+2 − 4𝑢𝑗+1 + 6𝑢𝑗 − 4𝑢𝑗−1 + 𝑢𝑗−2(17)
+ 𝑜(Δ𝑥 2 )Δ𝑥 4
次に，この微分値を直交格子において隣接セルの持つ値のみを用いて計算する方法を構築する．非構造格子において一般的に利用されている Green-Gauss 法もしくは Weighted-Least-Square(4)(WLSQ)法を用いると，Fig.2 のような等間隔格子上で勾配値は式(18)のように計算される．
𝑢𝑗+1 − 𝑢𝑗−1
(18)∇𝑢𝑗,𝐺 =
2Δ𝑥
ただし，添え字 G は Green-Gauss の方法で求めた値を表す．
また，同様にして隣接セルの値から 2 階の勾配値は以下のように与えられる．
𝑢𝑗+1 − 2𝑢𝑗 + 𝑢𝑗−1
(19)
∇2 𝑢𝑗,𝐺 =Δ𝑥 2
この値を用いると，式(20)～(22)のように勾配値が再構築できる．
なお，添え字 R は再構築された値を示す．
∇4 𝑢𝑗 =

式(9)を用いて 𝑢𝑗+1 を𝑛次精度で求めるためには，右辺をΔ𝑥 𝑛 の2項まで用いるとともに，𝑚階の勾配値を𝑛 − 𝑚次精度以上で求める必要がある．例えば，有限体積法における MUSCL でしばしば用いられてきた 2 次精度完全風上法での内挿では，右辺第 2 項まで用いたうえ，𝑢𝑗 − 𝑢𝑗−1
(10)
∇𝑢𝑗 =
+ 𝑜(Δ𝑥1 )Δ𝑥とする．このとき，式(9)に代入すると打切り誤差が𝑜(Δ𝑥 2 )となり
空間 2 次精度が達成できる．また，3 次精度風上バイアス法では右辺を第 3 項まで残し，

4
1 ∇𝑢(𝑗+1),𝐺 + ∇𝑢(𝑗−1),𝐺∇𝑢𝑗,𝑅 = ∇𝑢𝑗,𝐺 −
3
3
2

(20)

4
1 ∇𝑢(𝑗+1),𝐺 − ∇𝑢(𝑗−1),𝐺∇2 𝑢𝑗,𝑅 = ∇2 𝑢𝑗,𝐺 −
3
3
2Δ𝑥

(21)

∇2 𝑢(𝑗+1),𝐺 − 2∇2 𝑢𝑗,𝐺 + ∇2 𝑢(𝑗−1),𝐺
Δ𝑥 2

(22)

∇4 𝑢𝑗,𝑅 =2

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E04-3以上の手順によって，𝑗番目のセルにおける勾配値を，隣接セル(j-1，j+1 番目のセル)に格納されている変数の平均値および G-G
によって計算されたセル平均の勾配値のみを用いて高次精度で表すことが出来る．
ソルバーに実装する上では，計算負荷との兼ね合いから，式(9)で，(a) 4 次の勾配量を無視し 1 次の勾配量のみ高次精度化，(b) 4次の勾配量を無視し 1 次，2 次の勾配量を高次精度化，(c) 4 次の勾配量も実装し 1 次，2 次の勾配量を高次精度化，の 3 通りが考えられる．なお，演算量は少ない方から，(a)，(b)，(c)である．ここでは，内挿スキームの名称を，それぞれ 5 次精度風上バイアス
法(a)，(b)，(c)とする．
２．３．フーリエ解析
構築したスキームを格子幅 1 の等間隔格子上で 1 次元線型移流
方程式𝑢𝑡 = 𝑢𝑥 に適用し，フーリエ解析によって解像度を調べた．
手法の概略としては，まず，(23)𝑢(𝑥) = exp(𝑖𝑤𝑥)と置き，各スキームを用いて𝑢𝑥 の評価を行う．ただし，本手法は本来有限体積法による定式化であるが，これを 2 次精度の中心差
分と捉え，
𝜕𝑢
(24)
=𝑢 1−𝑢 1
𝑗+
𝑗−
𝜕𝑥
2
2
として差分評価を行い，
𝜕𝑢(25)= 𝑖𝑤 ∗ exp(𝑖𝑤𝑥)𝜕𝑥と整理する．
ここで，式(23)の解析的な微分は，
𝜕𝑢(26)= 𝑖𝑤 exp(𝑖𝑤𝑥)𝜕𝑥と表されるため，数値微分によって得られる修正波数𝑤 ∗ と本来の波数𝑤を比較し，勾配量の高次精度化による解像度への影響を調べた．なお，解析手法の詳細については文献(7)を参考されたい．
以上の手法に基づいて，2 次精度完全風上法，3 次精度風上バイ
アス法，および 5 次精度風上バイアス法(a)～(c)の各手法を用いて補間を行った場合について，1 階微分値の数値解の実部と虚部をFig.3, Fig.4 に示す．
図より，5 次精度風上バイアス法(a)～(c)を用いた場合は，3 次
精度風上バイアス法を用いた場合に比べ，より高い波数まで解像できることが分かる．また，5 次精度風上バイアス法(a)は (b)に比べ，虚部は劣るものの，実部では良い解像度を示している．5 次
精度風上バイアス法(c)は(b)に比べて実部，虚部ともに差がほとん
どなく，4 次の勾配量を計算する計算コストから考えると利点はないと言える．なお，以上の解析は線型問題に対して行ったものであり，実問題においては近似 Riemann ソルバーの影響が含まれる．

Fig.3 Real part of modified wave number

Fig.4 Imaginary part of modified wave number
２．４．不等間隔への対処
不等間隔部を含む格子に上記のスキームを適用したときに誤差がどのようになるかを考える．例えば，Fig.5 のようにハンギング
ノード近傍においてステンシルの先端が等間隔に存在しないときは，次のようにステンシルを考えていることになる．ここでは j-2のセルが等間隔位置に存在しないとすると，j-1 のセルに格納されている勾配値は最大 2 次精度なので，′𝑢𝑗−2= 𝑢𝑗−1 − ∇𝑢𝑗−1 Δ𝑥 + 𝑜(Δ𝑥 3 )(27)として等間隔位置に仮想の変数値を 3 次精度で設定し，その値を用いてステンシルを作成していることと同等となる．等間隔部分
ではこの仮想点での値と実際の j-2 番目のセルでの値が一致する．
また，セル j-1，j，j+1 の大きさが異なるときは，式(20)による勾配評価での誤差がどのようになるかを考える．ここでは各点が持つ勾配が正確な値であると仮定し，式(20)において右辺を点 j
周りでテーラー展開すると式(28)のようになる．

3

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E04-3

Fig.5 Imaginary stencil around hanging-nodes

Fig.6 Flux calculation on two-dimensional grid

4
∇𝑢𝑗,𝑅 = ∇𝑢𝑗
3

３．1 次元線型移流問題
３．１．問題設定
ここでは，構築したスキームの精度を検証するため，1 次元での密度波移流問題について調べる．計算領域は長さ2𝜋，両端は周期境界条件であり，この領域を 25～400 に等分割する．また，初期条件は計算領域内に 5 波長分の密度波が入るように，
𝜌 = 1.0 + 0.1 sin(5𝑥)(30)𝑢 = 1.0
𝑝 = 1.0とする．このとき，1 波長あたりのセル数は 5~80 となる．その他の計算条件は Tab.1 のようになっている．

Δ𝑥𝑗 + Δ𝑥𝑗+1
1− (∇𝑢𝑗 + ∇∇𝑢𝑗 () + 𝑜(Δ𝑥 2 ))6
2
(28)
Δ𝑥𝑗−1 + Δ𝑥𝑗
1+ (∇𝑢𝑗 − ∇∇𝑢𝑗 () + 𝑜(Δ𝑥 2 ))62
= ∇𝑢𝑗 + 𝑜(Δ𝑥1 )Green-Gauss，WLSQ 等の手法で勾配値を求めると最大で空間 2次精度となるが，式(28)を用いて再構築すると∇𝑢𝑗,𝑅 は空間 1 次精度に低下してしまう．そのため，式(20)をそのまま用いるのではなく，セルの大きさを用いて 1 次の誤差項を消すと以下のように
なる．
2Δ𝑥𝑗
2Δ𝑥𝑗
1
∇𝑢𝑗𝑅 = {1 + (+)} ∇𝑢𝑗
6 Δ𝑥𝑗 + Δ𝑥𝑗+1 Δ𝑥𝑗−1 + Δ𝑥𝑗
2Δ𝑥𝑗
2Δ𝑥𝑗
1
1
− (
) ∇𝑢𝑗+1 − () ∇𝑢𝑗−1(29)
6 Δ𝑥𝑗 + Δ𝑥𝑗+1
6 Δ𝑥𝑗−1 + Δ𝑥𝑗

Tab.1 Computational methods
Governing equation
Spatial discretization
Convection termTime integration
Gradient estimation= ∇𝑢𝑗 + 𝑜(Δ𝑥 2 )One-dimensional compressive Euler equation
Cell-centered finite volume method
SLAU
4th-order Runge-Kutta method
Green-Gauss method

なお，時間積分誤差の影響を排除するため，各計算において時間刻みを小さくしていき，解析解との誤差が一定値に収束することを確かめたうえで誤差評価を行う．

なお，等間隔部では式(29)は式(20)に一致する．
このように勾配を用いたステンシルの延長を行うことで，等間隔部分での精度向上を図りつつ，不等間隔部においても同じスキームで対応することができる．また，他形状格子を含むハイブリッド格子であっても，背景格子が直交等間隔格子であれば同手法を用いることができる．実際に G-G や WLSQ で求めた勾配をこの式に代入したときに誤差がどのようになるかは明確でないため，
この点については後に示す検証問題で調査する．

３．２．計算結果
1 波長あたりのセル数が 10 のときについて，各精度のスキームを用いて無次元時間 4 の間移流を行った結果を Fig.7～Fig.11 に示
す．

２．５．多次元問題への拡張
本手法を 2 次元以上の問題に適用する際には，式(1)における面積分の高次精度化が問題となる．ここで Fig.6 のような 2 次元セルにおいて x 方向の面への内挿を考えると，面中心の値𝑢𝐿 ，𝑢𝑅 は
上記の手法により空間 5 次精度で求まる．ここで，流束の面積分を正確に 4 次精度で行うとすると，面上でさらに面中心から複数点への内挿を行い，この点上でRiemann 流束を求める必要がある．
ここでは内挿のためにクロス微分項(𝜕 2 𝑞 ⁄𝜕𝑥𝜕𝑦など)も合わせ，3階までの 9 つ(3 次元では 16)の勾配値をセルに格納する必要があり，さらには流束評価の回数が 1 面あたり 2 回(3 次元では 4 回)に増える．面積分をさらに高次精度で行うと，演算量は飛躍的に増大する．
そこで本研究では補間手法のみを高次精度化した場合のスキー
ムの性質について調べた．
面積分の問題を回避する方法ついては，
今後有限差分法の考え方を応用するなどして取り組む予定である．

Fig.7 Density profile (2nd -order upwind)4

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E04-3
前述のように時間積分誤差は無視できるので，ここではΔ𝑥を小さくしたときの誤差の小さくなるスピードを見ることで，スキームの空間精度を確かめることができる．Fig.12 より， 5 次精度風上バイアス法の各手法を用いて1波長あたりのセル数10で計算した
結果は，3 次精度風上バイアス法で 1 波長あたり 20 の時よりも誤差が小さいことが分かる．さらに，5 次精度風上バイアス法(a)よりは(b)，(c)の精度が良くなっている．また，(b)と(c)では差がほとんどないことから，4 次の勾配量を計算して内挿式に組み込むことの影響は小さい．この傾向は 2.3 節で述べたフーリエ解析の結果とも一致する．また，実問題で考えると図中左側の，波長に対して格子数の少ない領域でどれだけ波を正確に捉えられるかが重要となり，右側の領域のようにもともと余裕のある部分での精度の持つ意味は少ない．以上の観点と，演算量の問題より以後の検証問題では 5 次精度風上バイアス法(a)，(b)のみについて調べる．

Fig.8 Density profile (3rd-order upwind-biased)

Tab.2 L1 norm of calculation error
2nd-order upwind

Cell number
25
50
100
200
400

L1
6.211E-02
3.709E-02
1.050E-02
2.637E-03
6.560E-04

accuracy
0.744
1.821
1.993
2.007

3rd-order upwind
biasedCell number
25
50
100
200
400

L1
6.132E-02
2.133E-02
3.193E-03
4.091E-04
5.137E-05

accuracy
1.523
2.740
2.964
2.993

5th -order upwind
biased (a)Cell number
25
50
100
200
400

L1
3.672E-02
2.495E-03
1.206E-04
6.757E-06
4.045E-07

accuracy
3.879
4.370
4.158
4.062

5th -order upwind
biased (b)Cell number
25
50
100
200
400

L1
3.720E-02
1.982E-03
6.647E-05
2.192E-06
6.875E-08

accuracy
4.231
4.898
4.922
4.995

5th -order upwind
biased (c)Cell number
25
50
100
200
400

L1
3.734E-02
1.956E-03
6.458E-05
2.088E-06
7.022E-08

accuracy
4.255
4.921
4.951
4.894

Fig.9 Density profile (5th-order upwind-biased (a))

Fig.10 Density profile (5th-order upwind-biased (b))

Fig.11 Density profile (5th-order upwind-biased (c))次に，計算結果と解析解のずれを L1 norm，つまり各点における差の絶対値の平均で表した結果を Tab.2 および Fig.12 に示す．
5

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E04-3

Fig.13 Initial value of density and computational grids
Fig.12 Spatial accuracy of each scheme
計算条件はTab.3に示した通りである．
空間補間法に関しては，
2 次精度完全風上法，3 次精度風上バイアス法および 5 次精度風上バイアス法(a)，(b)を用いる．5 次精度風上バイアス法(a)について
は 2.4 節で示した不等間隔への対処を行う場合と，行わない場合の 2 通りについて調べた．なお，時間刻み幅は無次元時間で 0.025
であり，600 [step]計算を行う．各手法においてさらに時間刻みを細かくしたときに結果がほぼ変化せず，ここでの時間積分誤差は十分小さいことを確認している．

４．2 次元等エントロピ渦移流問題
４．１． 問題設定
2 次元以上の問題における面積分の影響，およびハンギングノ
ードにおける計算の精度を調べるため，解析解のある 2 次元問題
について解析を行う．ここでは，等エントロピ渦を 45°方向に移流させ，解析解とのずれを調べる．一様流条件は(𝜌∞ , 𝑢∞ , 𝑣∞ , 𝑝∞ ) = (1, √2, √2, 1/𝛾)，比熱比𝛾 = 1.4である．
ここで，渦を含む流れ場の解析解は文献(8)を参考に，以下のよ
うに与える．
𝜌0 = [

1
(𝑇 −
𝑆∞ ∞

(𝛾 − 1)𝛤 2
8𝛾𝜋 2

2

e1−𝑟 )]

Tab. 3 Computational methods
Governing equation
Spatial discretization
Convection termTime integration
Gradient estimation1
𝛾−1

𝛤 1−𝑟 2(31)e 2 𝑦̅
2𝜋
𝛤 1−𝑟 2𝑣 0 = 𝑣∞ +e 2 𝑥̅
2𝜋
00𝑝 = 𝑆∞ (𝜌 )𝛾ただし，𝛤は渦の強さを表しており，ここでは 5 とする．𝑥̅ , 𝑦̅は渦𝑢0 = 𝑢∞ −Two-dimensional compressive Euler equation
Cell-centered finite-volume method
SLAU
4th-order Runge – Kutta method
Green-Gauss method

４．２．計算結果
無次元時間 5 ごとに計算結果を抽出し，一様流条件を差し引いたうえで足し合わせた結果を Fig.14～Fig.17 に示す．なお，渦は図中左から右方向に移流している．また，45°方向の断面を Fig.18に示す．

中心(𝑥𝑐 , 𝑦𝑐 )に対する座標を表しており，また𝑟 = √𝑥̅ 2 + 𝑦̅ 2 である．渦は一定速度で移流していくため，渦中心の座標は，(32)(𝑥𝑐 , 𝑦𝑐 ) = (𝑥𝑐0 , 𝑦𝑐0 ) + (√2, √2)𝑡と表される．さらに等エントロピ条件は𝑝/𝜌𝛾 = 𝑝∞ /𝜌∞ = 𝑆∞ として与えられる．
計算領域は一辺 32 の正方形とし，渦中心の初期位置は(𝑥𝑐0 = −12, 𝑦𝑐0 = −12)とする．格子幅は，𝑥 < −4かつ𝑦 < −4の領域のみ 0.25 とし，それ以外の部分が 0.5 となっている．以上の条件による密度の初期値，および計算格子を Fig.13 に示す．

6

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E04-3

Fig. 14 Density distribution (2nd-order upwind)

Fig.17 Density distribution (5th-order upwind biased(b),without treatment for non-uniform cells)

Fig.15 Density distribution (3rd-order upwind biased)

Fig.18 Cross section of density distribution
渦の移流結果を見ると，2 次精度完全風上法よりも 3 次精度風上バイアス法，3 次精度バイアス法よりも 5 次精度バイアス法(a)の方が低散逸で，より渦中心における密度低下の移流を正しく捉えられている．格子間隔の切り替え部は Fig.18 に示したように 2
つ目の山の内部に存在するが，その周辺で顕著な非物理的な散逸や反射，オーバーシュートが生じていないことが分かる．また断面を見ると，2 次精度完全風上法のみ大きく渦中心の位置のずれが生じている．これは Fig.3 に示した通り，2 次精度完全風上法では低周波域で実際よりも早く位相が進んでしまうことに起因する．
次に，計算結果と解析解の誤差の推移を調べる．ここでは各格子点における密度の誤差の絶対値を取り，足し合わせて全体の計
算誤差を L1 norm として評価した．時間積分による誤差は前述のように無視できるため，この誤差は各ステップにおける空間誤差の累積とみなせる．この結果を Fig.19 に示す．

Fig.16 Density distribution (5th-order upwind biased(a),without treatment for non-uniform cells)

7

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E04-3

Fig.21 Effect of treatment for non-uniform cells
Fig.19 Time history of L1 error
式(29)に示される対処を行った場合，不等間隔の界面を通過するときの誤差増加が抑えられる．
さらに，その後の誤差も減少し，最終値では約 6.7[%]小さくなっている．実際の流体計算を行う格子にはハンギングノードが多数存在するため，この対処は特に有効になると予想される．
最後に，計算時間について記す．本章で行った計算において，1次精度完全風上法(勾配計算なし)を用いて計算する時の計算時間を 1 とし，3 次精度風上バイアス法，5 次精度風上バイアス法(a)，(b)の各手法を用いたときの計算時間を Tab.4 に示す．各計算時間
は 5 回測定したときの平均値である．ただし，ここでは対流項の計算と時間積分にかかる時間のみを示している．
5 次精度風上バイアス法(a)の計算時間は 3 次精度風上バイアス法を用いたときの約 1.1 倍であり，3 次精度風上バイアス法との計
算時間の差は，1 次精度完全風上法と 3 次精度風上バイアス法の差，つまり，G-G によって勾配を求めるのに必要な時間よりも小さい．また，5 次精度風上バイアス法(a)と(b)で移流計算の結果があまり変わらないことを考えると，5 次精度風上バイアス法(a)は実際の問題に適用するうえで効率的であると考えられる．

図より補間法の精度が上がるにつれて，ステップごとの誤差の増加割合が小さくなることが分かる．さらに渦中心が格子間隔の切り替え位置に来たときを見ると，5 次精度風上バイアス法(a)，(b)ともに切り替え位置通過前後で若干の誤差増加が見られるが，全体としては 3 次精度風上バイアス法よりも良い計算結果が得られている．
次に，2.4 節で示した不等間隔への対処の有無について調べる．
ここでは，5 次精度風上バイアス(a)の手法について，式(29)に示される対処を行った場合の密度分布の推移を Fig.20 に示す．また，
この対処の有無による計算結果と解析解の誤差推移の比較を
Fig.21 に示す．

Tab.4 Averaged computation time
Scheme
1st-o upwind
3rd -o upwind biasedth
5 -o upwind biased (a)5th -o upwind biased (b)Computation Time
1
1.17
1.29
1.40

５．結論
非構造直交格子において高次精度でセル表面の変数値を補間する方法を構築した．非構造格子では隣接するセルの情報にしか簡
単にアクセスできないが，ここでは直交格子の利点を生かし，等
間隔領域において勾配値を高次精度化した．1 次元および 2 次元の検証問題において，新しく構築した手法を用いた場合，従来用
いられてきた 2 次精度や 3 次精度の手法を用いた場合よりも大幅に解像度が向上した．1 階勾配値のみを 4 次精度に向上させる方法は，補間スキームを完全な空間 5 次精度にする方法と比べ，検証問題の結果では大きな差異が見られなかった．また，ハンギン
グノードにおいては若干精度が低下するものの従来の手法より良thFig.20 Density distribution (5 -order upwind biased(a),with treatment for non-uniform cells))

8

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E04-3
い移流結果が得られ，また顕著な非物理的散逸や反射も生じないことが確認された．新しく提案した計算手法は，計算時間が従来の 3 次精度風上バイアス法と比べ 10[%]程度の増加に留まったことから，精度と演算量の両面から見てバランスのとれたスキームであると言える．

謝辞
本研究は，宇宙航空研究開発機構の委託研究「Ｈ25 年度非構造
格子における数値流束の高精度化」（JX-PSPC-380540）に基づき実施した．宇宙研究開発機構航空本部数値解析技術研究グループの橋本敦氏には論文作成にあたり貴重なコメントを頂いた．
また，基礎となる直交格子法については科学研究費助成事業，若手研究（B）（基金）
23760767 の成果を利用した．
ここに感謝の意を表す．

参考文献(1) 今村太郎，神園仁志，高橋悠一，“航空機空力性能評価を目的とした直交格子法ベース流体ソルバーの開発，”平成 23 年度航空宇宙空力班シンポジウム，2012．
(2) 神園仁志，今村太郎，“埋め込み境界法を用いた直交格子オイ
ラーソルバーの改良，”第 44 回流体力学講演会・航空宇宙数値シ
ミュレーション技術シンポジウム 2012．
(3) Shima, E., Kitamura, K. and Fujimoto K., “New gradient Calculation
Method for MUSCL Type Schemes in Arbitrary Polyhedra,” AIAA
2010-1081, 48th AIAA Aerospace Sciences Meeting Including the New
Horizons Forum and Aerospace Exposition, 2010.
(4) Mavriplis, D. J., “Revising the Least-square Procedure for Gradient
Reconstruction on Unstructured method,” AIAA 2003-3986, 16th AIAA
Computational Fluid Dynamics Conference, 2003.
(5) Cockburn, B. and Shu, C. W., “The Runge–Kutta discontinuous
Galerkin method for conservation laws V: multidimensional systems.,”
Journal of Computational Physics vol. 141 (1998): 199-224.
(6) Wang, Z. J., “Spectral (finite) volume method for conservation lawson unstructured grids. basic formulation: Basic formulation.,” Journal ofComputational Physics, vol. 178 (2002): 210-251.
(7) Lele, S. K., “Compact Finite Difference Schemes with Spectral-like
Resolution,” J. Computational Physics, vol. 103 (1992), 16-42.
(8) Corre, C., Hanss G. and Lerat, A.,“A residual-based compact scheme
for the unsteady compressible Navier-Stokes equations,” Computers &Fluids, vol. 34 (2005): 561-580.

9

Copyright © 2013 by JSFM