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

高次精度保存型 FR 法における一様流保持
Freestream preservation on a high-order conservative FR scheme阿部圭晃，東大院, 相模原市中央区由野台 3-1-1，E-mail：abe@ﬂab.isas.jaxa.jp

芳賀 臣紀，JAXA/JEDI, 相模原市中央区由野台 3-1-1，E-mail：haga.takanori@jaxa.jp
野々村 拓，JAXA/ISAS, 相模原市中央区由野台 3-1-1，E-mail：nonomura@ﬂab.isas.jaxa.jp
藤井 孝藏，JAXA/ISAS, 相模原市中央区由野台 3-1-1，E-mail：fujii@ﬂab.isas.jaxa.jp

Yoshiaki Abe, University of Tokyo, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：abe@ﬂab.isas.jaxa.jp
Takanori Haga, JAXA/JEDI, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：haga.takanori@jaxa.jp
Taku Nonomura, JAXA/ISAS, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：nonomura@ﬂab.isas.jaxa.jp
Kozo Fujii, JAXA/ISAS, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：fujii@ﬂab.isas.jaxa.jp
The evaluation of the symmetric conservative metrics are discussed in detail to satisfy a freestream preservation
using the conservative ﬂux reconstruction (FR) scheme on a three-dimensionally-curved stationary grid. In our
new evaluation, the conservative metrics(1)(2) is appropriately computed to avoid the freestream preservation
errors for any couple of the order of a solution polynomial and a shape function. The order of accuracy is also
validated for the computation of a convecting vortex on three-dimensional wavy grids. It is also numerically shownthat when the Radau polynomial is adopted for the modiﬁed function, the nonconservative metrics can preserve
the freestream althuogh the order of a solution polynomial is insuﬃcient. This is due to the cancellization of analiasing error in terms of Legendre polynomial when the ﬂux is reconstructed.

1.

はじめに

近年，不連続ガレルキン法（DG 法）に代表される非構造格子上で高い空間精度を有する手法の研究が進展し，複雑形状周りで効率の良い Large eddy simulation (LES) が実現しつつある．これらの非構造高次精度スキーム（DG，SD，SV，FR 法）は，空間を分割する計算要素（セル）内に高次精度化の為の自由度を持ち，セル境界の流束を近似リーマン解法等の数値流束で評価する．その為，総自由度を同じに保つ計算では，高次精度になれば計算コストが低くなる事が特徴である．さらに幾つかの問題では，従来の高次精度差分法と同程度か少ない自由度で同等の解像度を持つ例が報告されている(3) ．また，差分法を用いた複雑形状周りの圧縮性流れの数値計算では境界適合座標系を導入する事が多く，境界形状の精度は計算
格子の格子点数により直接的に決まる．一方，非構造高次精度スキームは各セルの境界形状を内部点と形状関数によって解析的に定義する為，事前に用意すべき計算格子の規模は著しく小さくなる（例えば 40 万点の格子での差分法と，5 万点の格子での境界形状 2 次精度の Flux
Reconstruction （FR）法では境界形状を表現する自由度が同じ）．この事は，今後エクサスケールでの高次精度差分法による流体計算で 100 億自由度の計算格子が必要となる際に，高次の形状関数を用いた FR 法では数千万点の計算格子で済む事を示唆しており，前処理の負担軽減の観点でも非常に有望なスキームと言える．本研究では特に，DG 法と同等の精度を有しより計算コストを抑えた FR 法(4) に着目し，特に保存型の支配方程式を解く
保存型 FR 法を扱う(5) ．
一般に，複雑形状周りの圧縮性流れの数値計算に境界適合座標系を導入すると，保存型の支配方程式を用いた
高次精度差分スキームでは一様流を保持しない．その為，

1

座標変換のメトリックを保存型の表式（保存型メトリック）に書き換え，これを解決する手法が提案されている(6)(7)(1)(8)(9)(2)(10)(11)(12)．保存型の支配方程式を用いたFR 法においても，セルの形状を定義する形状関数が高次である場合には同様の工夫が必要となる．Kopriva は SE
（Discontinuous spectral element）法において，差分法で導入された保存型メトリックを用いることで一様流保持を満たすスキームを提案した(13) ．彼は 3 次元の曲がった境界を持つ系で，解の次数が境界形状の次数の 2 倍未満
だと一様流を保持しない事を示し，保存型メトリックを用いると任意の次数で一様流を保持する事を解析・数値的に示している．しかし，以下の点は明確に示されていない．

• メトリック構築において保存量保存性の観点で注意
すべき点．
• 格子点の座標値から保存型メトリックを構築する詳
細な手順．
• 保存型メトリックを用いた際の数値解の精度検証．

図 1: Deﬁnition of grid, solution and ﬂux points in eachcell. The order of a shape function and solution polynoimial is N = 1 and M = 2, respectively.

Copyright c 2013 by JSFM

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

そこで本稿では，以下を新たにまとめる．

ここで，Q は保存量ベクトル，E, F, G はデカルト座標系
での流束関数である．
式 (2) を解く為に，各セルに内点（内部自由度）を導入し解を多項式近似する．一般座標の各方向毎に M + 1
個の Solution point （以下 SP：図 1 の緑色の点）を定義
し，各方向のラグランジュ多項式のテンソル積として保存量 Q̂ の第 n セル内の分布を与える．

1. 保存型 FR 法において一様流保持と保存量保存に必要な条件を示す（Sec. 3, Sec. 2.2）．
2. 一様流保持を満たすメトリックの評価方法を，保存量保存性及び補間関数の次数の関係に注意して詳細に説明（Sec. 4）．
3. 保存・非保存型メトリックを用いた数値解の精度検証を行い，所望の精度が得られるかを議論（Sec. 5）．

2.1

支配方程式の離散化

i,j,k=0

座標変換のヤコビアン J := |∂(x, y, z)/∂(ξ, η, ζ)| とメト
リック ξˆx = ξx J 等はこれらの微分を元に計算出来る．なお，図 1 は N = 1 の場合の 1 セルを示しており，GP は
青色のセル頂点にのみ存在する．N ≥ 2 の場合にはセル
境界上にも GP が定義されるが，実際に計算格子を準備する際には頂点の GP のみを定義すればよい（境界上のGP は計算中に定義されるものであり，計算空間で等間隔に置けば良い）．以下では簡単の為，事前に定義する
頂点の GP と，計算中に定義されるセル境界上 GP を同じ GP で表す．計算空間における圧縮性 Euler 方程式は次のように書ける1 ．

Q̂ = QJ,
1

Ê = ξˆx E + ξˆy F + ξˆz G, . . .

SP
φMn;p (ξ)Ên;p,q,r ,(5)∂ξ ÊnSP (ξ) = ∂ξ I M [ÊnSP ](ξ) ≡

M
∑

SP
dφM
n;p (ξ)Ên;p,q,r , (6)p=0φMn;p は以下で決まる M 次ラグランジュ多項式であり，
M
dφMn;p = dφn;p /dξ とする．
φMn;p (ξ) =M∏
m=0,m6=p

ξ − ξm.
ξp − ξm(7)

この流束 ÊnSP は隣接セルの情報を反映せず，セル境界で不連続な分布となる．そこで，セル境界で共通の値を採るような分布 ÊnC を再構築する．まず，セル境界に各方向の SP を通る直線とセルを囲む面の交点に Flux point（以下 FP：図 1 の赤色の点）を置く．式 (4) を用いて境界の両側のセルから FP における解 Q̂L , Q̂R を外挿し，近似リーマン解法を用いて Ê com を風上的に求める．そして，FP で Ê com を通るように元の M 次の ÊnSP を修正
し ÊnC を得る.
[]
com
ÊnC (ξ) = ÊnSP (ξ) + Ên−1/2− ÊnSP (−1) gL (ξ)[]
com+ Ên+1/2− ÊnSP (+1) gR (ξ) (8)

N
N
N
Mn;i
(ξ)Mn;j(η)Mn;k(ζ)rn;i,j,k (1)∂τ Q̂ + ∂ξ Ê + ∂η F̂ + ∂ζ Ĝ = 0,

M
∑
p=0

rn (ξ, η, ζ) = I N [rn ]N∑

p,q,r=0ÊnSP (ξ) = I M [ÊnSP ](ξ) ≡計算領域は六面体セルで分割されているものとする．
デカルト座標系 (x, y, z) で表される物理空間の全てのセルを，それぞれ一般座標 (ξ, η, ζ) で表される計算空間の標準セル Es = {ξ, η, ζ| − 1 ≤ ξ, η, ζ ≤ 1} に変換する．
物理空間の第 n セルの形状を表す為に，一般座標の各方
向毎に N + 1 個ずつの Grid point (以下 GP) を定義し，それらの座標 rn;i,j,k に対応する各方向 N 次の形状関数
MnN を用いる（図 1）．

=

M
M
φM
n;p (ξ)φn;q (η)φn;r (ζ)Q̂n;p,q,r

本研究では SP にガウス点を用いる．M 次の近似多項式（各方向毎に M + 1 個の SP）を用いると，空間 M + 1
次精度のスキームとなる．
各 SP における流束の発散を計算する手順を説明する．
一般性を失わずに Ê の ξ 方向微分のみを取り上げる（η =ηq , ζ = ζr と固定）．第 n セル内の分布 ÊnSP , ∂ξ ÊnSP はSP
式 (4) に倣い，SP での離散値 Ên;p,q,rを用いて M 次多項式で表現出来る．

4. 従来の非保存型メトリックで修正関数に Radau 関数を用いると，解の次数が形状関数の次数の 2 倍未満
でも一様流を保持する事を，ルジャンドル多項式によって解析的に説明（Sec. 6）．
Flux Reconstruction 法

M
∑

(4)

項目 3 の数値計算の結果，従来の非保存型のメトリック
でも修正関数に Radau 関数を用いると，解の次数が形状関数の次数の 2 倍未満でも一様流を保持するという新しい結果を得た．そこで，以下を本稿の内容に追加する．

2.

Q̂n = I [Q̂n ] ≡M

ここで，左右のセル境界は n ∓ 1/2 で示す．gL (ξ) は左端（ξ = −1）で 1，右端（ξ = +1）で 0 となる M + 1 次多項式であり，修正関数(4) と呼ばれる．本稿では以下で定義される関数を採用し，gL = gGa,M +1 とする．

(2)

M
M −1
RR,M +
RR,M −1 , (9)2M − 1
2M − 1(−1)M
RR,M = gDG,M =(PM − PM −1 ),(10)2gGa,M =

(3)

以下では，物理量 ψ の x 方向の偏微分表記として ∂x ψ 及び ψx を併用する．

2

Copyright c 2013 by JSFM

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

PM は M 次のルジャンドル多項式である．また，gR は
gL の原点対称な関数（gR (ξ) = −gL (−ξ)）で与えられる．RR,M は右ラダウ多項式と呼ばれ，修正関数 gL =gDG,M +1 = RR,M +1 として用いる事も出来る（Sec.6. で使用）．この場合，DG 法に等しいスキームとなる事が知られている(4) ．
2.2 保存量保存性について
保存型の支配方程式を用いるメリットの 1 つに，計算領域内での保存量積分値が一定に保たれるという保存量
保存性がある(2) ．FR 法では離散化が適切でないとこの性質が失われる場合がある．
まず，各セル内での保存性（局所保存性）を議論する．
FR の離散化は各座標方向に独立して行われる為，1 次元で考える．第 n セル内の保存量が時間発展に対して保存するには次式を満たす必要がある．
∫
com
com
∂ξ ÊnC dξ = Ên+1/2− Ên−1/2
.
(11)
Ωn

式 (8) で得られる流束 ÊnC は境界を含むセル内で連続な
M 次の多項式であり，式 (11) は満足する．SP としてガ
ウス点を用いている為，式 (8) の ÊnSP (±1) は式 (6) のÊnSP により外挿された値である．この時，保存量 Q を外挿して境界の流束を計算すると内部の近似多項式との不整合を生じ，局所保存性が失われる事に注意する．
次に，計算領域内での保存性（保存量保存性）を議論する．任意の閉領域 ∀ Ω における保存量保存性が満たされるには次式を満たす必要がある．
∑ ∫∂ξ ÊnC dξ = Ê∂com− Ê∂com
.
∀ Ω;n
∀ Ω;n
max +1/2
min −1/2n∈∀ Ω

FR 法における一様流保持の条件
一様流保持は差分法でしばしば問題となるが，FR 法や DG 法でも任意格子形状，特に高次の形状関数で表現されたセルに対する一様流保持(13) の詳細な議論は十分でない．本節ではまず，FR 法において一様流保持が満たされる為に必要な条件を整理する．

3.

3.1 一様流保持と面積保存則
一様流計算時の再構築流束 Ê C (ξ) を書き下す．なお以
下は全て第 n セル内の議論であり，下付き添字の n は一部省略する．
[]com
Ê C (ξ) = Ê SP (ξ) + Ên−1/2− Ê SP (−1) gL (ξ)[]
com+ Ên+1/2− Ê SP (+1) gR (ξ)= E uni ξˆxC (ξ) + F uni ξˆyC (ξ) + Guni ξˆzC (ξ)

Ωn

(14)

ここで，E uni , F uni , Guni は x, y, z 方向の一様流流束，
ξˆxC = {ξx J}C 等はメトリックである．一様流保持を満たすには，SP において ξˆxC 等が以下の面積保存則（: SCL）
を満たす必要がある．

(12)
右辺は閉領域 ∀ Ω の境界積分に相当する．各セルで局所保存性が成り立つ場合，下記のように書き換えられる．
∑
com
com
)
(Ên+1/2− Ên−1/2n∈∀ Ω= Ê∂com− Ê∂com
.
∀ Ω;n
∀ Ω;n
max +1/2
min −1/2

構築すれば良い．このように隣り合うセルで GP を共有させる為に，各方向に分布させる GP の端点は必ずセルの境界に位置していなければならない．以下で新たに提案する一様流保持が可能なメトリックはこれらの性質を満たすように構築しており，保存量保存性が保証されている．
なお，FP において共通のメトリックを使用せずに一様流を保持するスキームを構築する事は可能であるが，その場合一般には保存量保存性が満たされない事を強調しておきたい．各セル内で独立に（局所的に）議論される
一様流保持と，隣り合うセルの持つ値まで考慮して議論される保存量保存性の概念を両立する為には特別な工夫が必要である．

(13)

com
もし各セル境界（FP）の Ên±1/2について，隣り合うセル同士で同じ値を共有すれば式 (13) 左辺の総和は閉領域∀Ω 境界の値のみが残る．ゆえに，各セル境界（FP）でcom
Ên±1/2の値が一意に定まる事は保存量保存性が成り立つ
為の十分条件である．これは直交等間隔格子では自明に
com成り立つ．しかし格子形状が曲線の場合には，Ên±1/2に含まれるメトリックもセル同士で共有する FP で同じ値を持たなければならず自明でない．本研究では，差分法における特性インターフェイス条件で提案されたメトリッ
クの評価方法の考え(14)(15) を拡張し，FP でのメトリックを構築する際に隣り合うセル同士で同じ離散値から補間して与える事を提案する．例えば，ξx J = yη zζ − yζ zη を
n + 1/2 番目の FP で定義する際には，n, n + 1 番目の両セルで共有する GP の値（y, z ）のみを用いて補間関数を∂ξ ξˆxC + ∂η η̂xC + ∂ζ ζ̂xC = 0,(15)∂ξ ξˆyC + ∂η η̂yC + ∂ζ ζ̂yC = 0,(16)∂ξ ξˆzC + ∂η η̂zC + ∂ζ ζ̂zC = 0,(17)

以下では x 方向の面積保存則 (15) に注目する．

3.2 面積保存則の成立条件
先に示したメトリック ξˆxC を以下では再構築メトリッ
クと呼ぶ事にする．式 (14) に従うと，ξˆxC は[]
ξˆxC (ξ) = ξˆxSP (ξ) + ξˆxGP (−1) − ξˆxSP (−1) gL (ξ)[]+ ξˆxGP (+1) − ξˆxSP (+1) gR (ξ) (18)と書ける．ここで，ξˆxSP は SP でのメトリックの離散値
から解の次数で構築した補間関数，ξˆxGP は GP でのメト
リックの離散値から構築した補間関数である．引数の ±1は，各セルの左右の FP を示す．注意すべきは，FP におけるメトリックを ξˆxGP (±1) とし，隣り合うセル同士で共3

Copyright c 2013 by JSFM

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

細に述べる．i, j, k, l, m, n は GP，p, q, r は SP（：ガウス点) を指す．まず，通常用いられる NC を用いたメトリッ
ク ξˆxC は M ≥ 2N の時のみ面積保存則を満足する．次に，SC を導入すると M ≥ N の時に面積保存則を満足する．
さらに評価を工夫した SCHGP を用いると任意の M, Nのセットに対して面積保存則を満たす．
4.1 非保存型メトリック（NC: ξˆx = yη zζ − yζ zη ）有する FP で同じ値を採るようにしなければ保存量保存性が破綻する事である．式 (18) より面積保存則 (15) を書き下すと以下の通りである．
式 (15) 左辺 ⇐⇒ ∂ξ ξˆxSP + ∂η η̂xSP + ∂ζ ζ̂xSP
SP part[]+ ξˆxGP (−1) − ξˆxSP (−1) dgL (ξ)/dξ + . . .
Modified part

通常，メトリックには NC が用いられる．

(19)

式 (19) の下線部（SP 及び Modiﬁed part）が SP におい
Step 1: N 次の補間関数 y 等 y, z の補間関数を，格子て 0 となる条件を次にまとめる．
点（GP）の離散値に基づき N 次の形状関数で定義する．
メトリックの計算時に以下のいずれかの条件を満たせ
N∑
ば SP part は SP において 0 となる．
N
I
[y]
=
MiN (ξ)MjN (η)MkN (ζ)yi,j,k .
(23)
SP parti,j,k=0(A-1) 解の次数 M が形状関数の次数 N の 2 倍以上

I N [z] =（2N ≤ M ）N∑

MiN (ξ)MjN (η)MkN (ζ)zi,j,k .

(24)

i,j,k=0(A-2) 対称保存型メトリックを用いる（補遺 B に示
Step 2: N 次の補間関数 yη 等 yη , yζ , zη , zζ の補間関数
す次数の制限有り）（N 次）を，格子点の座標値と形状関数を用いて計算する．

N
条件（A-1）が成立すると ∂ξ ξˆxSP (ξ) 等は数値計算上，∑
N∂η I [y] =MiN (ξ)dMjN (η)MkN (ζ)yi,j,k , (25)
Chain rule に基づいた解析的な変形の結果と等しくなる
i,j,k=0
為，SP で恒等的に式 (15) が成立する．しかし，解の次数N∑
M が形状関数の次数 N の 2 倍未満だと非保存型メトリッ∂ζ I N [y] =MiN (ξ)MjN (η)dMkN (ζ)yi,j,k , (26)クは数値的に式 (15) を満たさない．一方，条件（A-2）のi,j,k=0対称保存型メトリック（: SC）は以下のような表式であN∑る（TA , TB , TC は下式で定義する）．
∂η I N [z] =M N (ξ)dM N (η)M N (ζ)zi,j,k , (27)iξˆx = {(yη z − zη y)ζ − (yζ z − zζ y)η }/2 = TAζ − TBη (20)∂ζ I N [z] =η̂x = {(yζ z − zζ y)ξ − (yξ z − zξ y)ζ }/2 = TBξ − TCζ (21)

これらは解析的には通常のメトリック表現（非保存型メ
トリック: NC）と等価だが，離散化・補間を伴うと一般に異なる．対称保存型メトリックを用いると補遺 B の条
件下で，微分と補間の可換性から SP part は 0 となる．
なお，メトリック計算過程の次数の変化は次節で詳細に議論する．
メトリックの計算時に以下の条件を満たせば Modiﬁed
part は 0 となる．

Modiﬁed part=N∑

MiN (ξ)MjN (η)dMkN (ζ)zi,j,k ,(28)

Step 3: 2N 次の補間関数 ξˆxGP ;2N
ξˆx の補間関数（2N次）を，格子点での yη 等から計算する．
ξˆxGP ;2N = ∂η I N [y]∂ζ I N [z] − ∂ζ I N [y]∂η I N [z](29)

これは N 次の形状関数で定義される格子形状の解析的なメトリックに一致する為，以下では ξˆxGP ;2N を真のメトリックと呼ぶ．

Step 4: M 次の補間関数 ξˆxSP ;Mξˆx の 補 間 関 数（MGP;2N次）を，先の ξˆ
の SP における離散値から再構成
x

する．



なお，先の SP part における条件（A-1）が成立する時には，本条件（B）も同時に成立する．
4.

ki,j,k=0ζ̂x = {(yξ z − zξ y)η − (yη z − zη y)ξ }/2 = TCη − TAξ (22)(B) ξˆxGP (±1) = ξˆxSP (±1) か つ η̂xGP (±1)
η̂xSP (±1) かつ ζ̂xGP (±1) = ζ̂xSP (±1)

ji,j,k=0ξˆxSP ;M = I SP ;M [ξˆxGP ;2N ]=メトリックの計算方法

M
∑

M
M
ˆGP ;2N
φM
p (ξ)φq (η)φr (ζ)ξx;p,q,r ,(30)

p,q,r=0

ここでは，非保存型（NC），対称保存型（SC），対称保存型高次メトリック（SCHGP）の 3 つを導入し，SPで用いる ξˆxSP (ξ) と FP で用いる ξˆxGP (ξ) の計算方法を詳以上で ξˆxGP ;2N と ξˆxSP ;M が求まり，ξˆxC の式 (14) を計算出
来る．M < 2N の場合，ξˆxSP ;M は真のメトリック ξˆxGP ;2N

4

Copyright c 2013 by JSFM

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

る TA , TB の微分と，流束の微分（面積保存則 (15) での微分）の可換性により面積保存則を満たす事が出来，2N 次の真のメトリックを構築する必要がない（補遺 C 参照）．
これにより，NC だと M > 2N の時にのみ面積保存則を満たすという条件を，M > N にまで緩める事が出来る．
Step 1,2 までは NC と同じである．

TA , T B のStep 3: 2N 次の補間関数 TAGP ;2N , TBGP ;2N
補間関数（2N 次）を，格子点の座標値と形状関数を用いて計算する．
2TAGP ;2N = ∂η I N [y]I N [z] − ∂η I N [z]I N [y](31)

2TBGP ;2N = ∂ζ I N [y]I N [z] − ∂ζ I N [z]I N [y].

(32)

Step 4: N 次の補間関数 ξˆxGP ;N ξˆx = TAζ − TBη の補間関数を N 次の形状関数を用いて定義する．
ξˆxGP ;N = ∂ζ I N [TAGP ;2N ] − ∂η I N [TBGP ;2N ]N{∑
GP ;2NMiN (ξ)MjN (η)dMkN (ζ)TA;i,j,k=

i,j,k=0}
GP ;2N− MiN (ξ)dMjN (η)MkN (ζ)TB;i,j,k. (33)

図 2: Costruction of NC in the case with N = 1 andM = 2.
（2N 次）を低次（M 次）で再構成する（Step 4）ので誤差を生じる2 ．この時 ξˆxSP ;M は面積保存則を満足せず（SPpart6= 0），かつ FP で ξˆxGP ;2N 6= ξˆxSP ;M となる（Modiﬁed
part6= 0）．ゆえに，NC を用いると M ≥ 2N の時に限り面積保存則が満足される（表 1）．これにより，例えば形状関数が 6 次の時には解の次数を 12 次にしなければな
らない等，形状関数の次数が高い程飛躍的に実装が難し
くなる．また，移動変形格子を使用する際にはメトリッ
クの計算に 3 次の非線形項（例えば xt yζ zη ）が現れる為，M ≥ 2N + Nt （N は空間方向の次数，Nt は時間方向の次数）というさらに厳しい要件となる．

Step5: M 次の補間関数 ξˆxSP ;M ξˆx の補間関数（M 次）を，先の ξˆGP ;N の SP における離散値から再構成する．

表 1: Condition for a freestream preservation using NC(◦ and • indicate zero and nonzero, respectively)

Cases

SP partModiﬁed partM < 2NM ≥ 2N•◦ [ by (A-1) ]•◦ξˆxGP ;N では 2N 次の TA , TB の離散値を N 次で補間し微分する為，NC で与えた真のメトリック ξˆxGP ;2N には必ず一致しない．
なお，ξˆxGP ;N に現れる微分は隣り合うセルで共有する境界面内の微分（η, ζ 方向）である．また被微分項である TA , TB の離散値には，隣り合うセルで共有する GP において定義された値のみを用いる．これらの事から，隣り合うセルで共有する境界面内で定義される ξˆxGP ;N は一通りに定まり，FP におけるメトリックの値は隣り合うセルで共通の値となる．これは Sec.2.2 で議論した保存量保存性を保証する為の十分条件である．また隣り合うセル
で GP を共有する為に，各方向に分布させる GP の端点は必ずセルの境界に置かなければならない．

x

ξˆxSP ;M = I M [ξˆxGP ;N ]
=

M
∑

M
M
ˆGP ;NφMp (ξ)φq (η)φr (ζ)ξx;p,q,r ,(34)

p,q,r=0

4.2

M < N の場合，ξˆxSP ;M は ξˆxGP ;N（N 次）を低次（M 次）
で再構成するので誤差を生じる3 ．この時，補遺 C より
ξˆxSP ;M は面積保存則を満たさず（SP part6= 0），かつ FPで ξˆxSP ;M 6= ξˆxGP ;N となる（Modiﬁed part6= 0）．SC が面積保存則を満たす条件は以下のようになる（表 2）．

対称保存型メトリック（SC: ξˆx = {(yη z − zη y)ζ −(yζ z − zζ y)η }/2)
ここでは，対称保存型メトリック（SC）の計算方法を述べる．SC を用いると，メトリックを構築する際に用い
2
SP ;M
GP ;2N
SP では ξ̂x= ξ̂xだが，他の場所では一致しない（補遺 B 参
SP ;M
GP ;2N照）．ゆえに ∂ξ ξ̂x
6= ∂ξ ξ̂xであり，面積保存則の各項が解析的なそれと一致しない事を意味する．

3
GP ;N
SP ;M
M < N の時には SP 以外では ξ̂x
6= ξ̂xとなる為，式 (33), あるいは補遺 C の式 (63) の形には書けない．

5

Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
C04-2を考える．Step 1,2 までは NC と同じ，Step 3 は SC と同じである．

表 2: Condition for a freestream preservation using SC(◦ and • indicate zero and non-zero, respectively)

Cases

SP partModiﬁed part

M <N
M ≥N

•
◦ [ by (A-2) ]•◦本手法では解の次数 M が十分に高くとも（M ≥ 2N ），ξˆxSP ;M = ξˆxGP ;N となるだけで，NC で与えた真のメトリッ
ク ξˆxGP ;2N は得られない．ゆえに，解の次数 M を高くしてもメトリックの次数は N 次以上にならず，得られる
数値解はメトリックの次数（N ）以上の精度にならない（Sec. 5.2）．このように一様流を保持する保存型メトリックを用いても，メトリックの次数まで考慮に入れ適切に評価しないと解の精度が出ない場合がある．
図 4: Costruction of SCHGP in the case with N = 1
and M = 2.

Step 4: M 次の補間関数 ξˆxGP ;Mξˆx = TAζ −TBη を M
次の補間関数を用いて定義する．ここで，hi, hj, hk は
(M + 1) × (M + 1) × (M + 1) 個の定義点（以下 HGP と呼ぶ）を表し，SP と一致する必要は無い．しかし保存量
保存性を考慮し，各方向の端点はセル境界に置かなければならない4 （図 4 の三角形の位置を参照）．ゆえに，各セルの頂点の GP には必ず HGP が置かれる．
ξˆxGP ;M = ∂ζ I M [TAGP ;2N ] − ∂η I M [TBGP ;2N ]
=

M
∑

{
GP ;2NM
M
M
Mhi(ξ)Mhj(η)dMhk(ζ)TA;hi,hj,hkhi,hj,hk=0
GP ;2NM
MM− Mhi(ξ)dMhj(η)Mhk(ζ)TB;hi,hj,hk}(35)

ξˆxGP ;M では 2N 次の TA , TB の離散値を M 次で補間し微分する為，M ≥ 2N の時には NC で与えた真のメトリッ
ク ξˆxGP ;2N と一致する．

図 3: Costruction of SC in the case with N = 1 andM = 2.

4.3

対 称 保 存 型 高 次 メ ト リック（SCHGP: ξˆxStep 5: M 次の補間関数 ξˆxSP ;Mξˆx = TAζ − TBη を M
次の補間関数を用いて定義する．

=
ξˆxSP ;M = I M [ξˆxGP ;M ]{(yη z − zη y)ζ − (yζ z − zζ y)η }/2）
先の SC では M < N の場合に面積保存則が満たされない．これは GP で構築したメトリック ξˆxGP ;N を解の次
数 M で補間再構築すると誤差を生じる事による．また，解の次数 M を十分高くしてもメトリックの次数は N 次
以上にはならず，数値解の精度に影響が出る（Sec. 5.2）
そこで，GP で構築するメトリックを解の次数 M に合わせ，任意の N, M に対して面積保存則を満たす評価手法

6

=

M
∑

M
M
ˆGP ;MφM
p (ξ)φq (η)φr (ζ)ξx;p,q,r ,(36)

p,q,r=0ξˆxSP ;M では ξˆxGP ;M （M 次）を M 次で補間する為，常にξˆxSP ;M = ξˆxGP ;M である．補遺 C より ξˆxSP ;M は面積保存則
4
FP において，隣り合うセル同士で同じメトリックを共有する為である（詳細は SC の Step4 及び Sec.2.2 参照）．

Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
C04-2を満足し（SP part は 0），かつ FP でも ξˆxSP ;M = ξˆxGP ;Mとなる（Modiﬁed part は 0）．SCHGP が面積保存則を満たす条件は以下のようになる（表 3）．
zy表 3: Condition for a freestream preservation using SC(◦ and • indicate zero and non-zero, respectively)

Cases

SP partModiﬁed partM < 2NM ≥ 2N◦ [ by (A-2) ]◦ [ by (A-1) ]◦◦図 5: Computational grid and initial ﬂow (convecting
vortex)5.1 一様流保持性
初期解 u0 = (u0 , v0 , w0 ), p0 , ρ0 として，x 方向の一様流を以下のように与える．

以上より，各メトリックを用いた場合の一様流保持は解の次数 M ，形状関数の次数 N に対して表 4 の通りまとめられる。

u0 = 0.10,

Order of SP M

NC

SC

SCHGP

M <N
N ≤ M < 2N

•
•

•
◦

◦
◦

2N ≤ M◦◦◦v0 = w0 = 0.0,

p0 /γ = ρ0 = 1.0,(40)

格子点数は各方向に Nmax = 10（自由度は解の次数によって異なる）とした．形状関数の精度が N = 2 の時，50 ステップ計算時の流れ場の L∞ ノルム（v, w の最大誤差）を表 5 に示す．これにより，一様流を保持するメト

表 4: Expected characteristics of a freestream preservation using each metric evaluation (◦ satisﬁes afreestream preservation; • does not satisfy a freestream
preservation)

5.

x

表 5:
L∞ norm of a freestream error(max{|v − v0 |/u0 , |w − w0 |/u0 }) using each metricevaluation with gGa modiﬁed fuction. The order of ashape fuction is N = 2.
Order of SP M
1 (M < N )2 (N ≤ M < 2N )
3 (N ≤ M < 2N )
4 (2N ≤ M )検証計算

前節で定義した NC,SC,SCHGP を用いて，一様流保持と移流渦の検証計算を行う．計算対象は圧縮性非粘性
流れとし，支配方程式には音速と一様流密度で無次元化された式 (2) を用いる．気体は理想気体とし，比熱比は
γ = 1.4 とする．本計算には，−5 ≤ x, y, z ≤ 5 なる立方体内の直交格子 x0 , y0 , z0 を式 (37)-(39) で 3 次元的に歪めた計算格子を用意し，その１つを図 5 に示す．n =4, A = 0.4, L0 = 10 とし，∆x0 = ∆y0 = ∆z0 はセル頂点（Node point）間の格子幅とする．流れ場には周期性を仮定し，計算領域の境界では各方向に周期境界条件を課した．

NC
5.2 × 10−2
1.9 × 10−2
4.4 × 10−3
1.5 × 10−13SC
4.9 × 10−2
6.8 × 10−14
1.4 × 10−13
2.8 × 10−13

SCHGP
1.3 × 10−14
6.8 × 10−14
1.5 × 10−13
2.8 × 10−13

リックと次数の関係は前節の表 4 の通りである事が数値的にも確かめられた．

5.2 移流渦
次に，2 次元渦(16) を移流させた計算結果を以下に示す．
移流渦は以下のように定義した．
u = U∞ − C[]nπ(k − 1)∆y0nπ(l − 1)∆z0x = x0 + ∆x0 A sin
sin
,
L0
L0
(37)
[
]
nπ(j − 1)∆x0nπ(l − 1)∆z0y = y0 + ∆y0 A sin
sin
,
L0
L0
(38)
[
]
nπ(j − 1)∆x0nπ(k − 1)∆y0z = z0 + ∆z0 A sin
sin
.
L0
L0
(39)

v=C

ρC 2
(y − yc )−r2, p = p∞ −exp
exp(−r2 ),
2
Rc
2
2Rc2
(41)

(x − xc )[(x − xc )2 + (y − yc )2 ]
−r2 2
,r =
exp
2
Rc
2
Rc2
(42)

C = U∞ ,(43)

なお，Rc = L0 /10,  = 0.02, U∞ = 0.5, xc = yc =0, ρ = 1.0 とした．移流した渦が一周して元の位置に戻ってくるまで計算し，z 一定の面の圧力分布を図 6 に示す．
図 6 と一様流保持の表 5 を比較すると，一様流保持を満たさないケースでは渦を全く保持出来ていない．一方，一7

Copyright c 2013 by JSFM

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

様流保持を満たすケースでは渦を保持出来ている．特に解の次数が M = 4 の時，厳密なメトリックを表現する
NC と SCHGP は，そうでない SC（Sec. 4.2）よりも渦を良く保持出来ている．これは，一様流保持を満たしていても，メトリックの次数が解の次数よりも低いと渦を良く保持出来ない事を示している．

表 6: The number of node points prepared for the convergence study on each accuracy M.

NDOF \ SP1

3

2

O(30 ): Coarse

3

3

15

10

O(403 ): Medium
O(503 ): Fine

183
243

123
163

3

4

3

7

63

93
123

73
103

次に，Swirl velocity: v の L2 ノルムを図 7 に載せる．
ここでは NDOF = O(303 ), O(403 ), O(503 ) の格子を用意し，自由度に対する収束性を見た．なお，表 6 には各次数での計算に用いた格子点数（セル頂点の数：Node point）を示した．

(a) NC (M = 1)(b) SC (M = 1)(a) NC(c) SCHGP(M = 1)(b) SC(d) NC (M = 2)(e) SC (M = 2)(f) SCHGP(M = 2)(g) NC (M = 3)(h) SC (M = 3)(i) SCHGP(M = 3)(j) NC (M = 4)(k) SC (M = 4)(l) SCHGP(M = 4)(c) SCHGP

図 7: The order of accuracy using NC, SC and SCHGP
for metric evaluation.

図 6: The contours of pressure distribution on vertical(z-constant) plane are shown.

一様流保持を満たさないケース（図 7(a)：NC の p1-p3,
図 7(b)：SC の p1）は解の収束性が悪い．一方，全ての次

8Copyright c 2013 by JSFM

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

数で一様流保持が満たされる SCHGP では，解の次数が高くなるにつれて収束性が良くなる．ここで，一様流保持が満たされても解の収束性が悪いケース（図 7(b)：SCの p3, p4）に注目する．Sec. 4.2 で触れたように SC で表現されるメトリックは，解の次数が M > 2 の場合でも形状関数の次数（N = 2）以上の次数を持たない．ゆえに，これらのケースでは数値解の精度がメトリックの精度で制限されていると考えられる（p3,p4 においても p2と同じ収束性しか持たない）．このように，保存型メト
リックにより一様流を保持するスキームを用いても，メ
トリックの次数が解の次数よりも低い場合には数値解の精度が出ない場合がある事が分かった．

Sec.4. 及び 5. で説明したように gGa を用いた NC がM < 2N の時に一様流を保持しない理由は，SP の次数で構築されるメトリックの補間関数 ξ SP （M 次）が格子形状の解析的な ξˆx（2N 次）を表現出来ない事にある．この為，SP でのメトリック ξˆxSP は面積保存則を満足せず（SP part6= 0），かつ FP において ξˆxGP 6= ξˆxSP （Modiﬁed
part6= 0）となり，再構築メトリック ξˆxC は面積保存則を満たさない．しかし SP part と Modiﬁed part のそれぞれが 0 でなくとも，互いに打ち消し合って 0 となれば再構築メトリック ξˆxC は面積保存則を満たす．修正関数に
gDG を用いた際，解の次数 M = 2N − 1, 2N − 2 の時には，このような特殊な状況が生じる．

なお，いずれも超収束の傾向にあるが，このような傾向は FR 法や DG 法で線形問題を解く際に一般に見られるものである．今回の検証計算では移流渦を取り扱ったが，何らかの原因でこの傾向が生じたと考えられる．
最後に，諸言で指摘した自由度に対する計算時間について触れる．表 7 は，自由度 303 における SCHGP の各次数での計算時間（計算機には Intel(R)Xeon5450(3.0GHz) を使用，自動並列における ifort の最適化オプションを O3）である．表 7 から，同じ自由度での計算時間は次数が高い程短くなっている事が分かる．

説明を簡単化する為，1 次元関数（ξ のみの関数）を仮定する5 ．まず，Sec.4. の NC における Step 3 の ξˆxGP ;2Nを想定し，2N 次の関数 f GP ;2N (ξ) を考える．さらに Step
4 の ξˆxSP ;M を想定し，M + 1 個の Gauss 点におけるf GP ;M (ξ) の離散値から，f GP ;M (ξ) の M 次補間多項式を構築する（M < 2N とする）．このように得た補間多項式 f GP ;2N ,f SP ;M は，ルジャンドル多項式 Pn (ξ) で展開して以下のように書ける．

f GP ;2N (ξ) =2N∑

Pi (ξ)bi ,(44)

Pi (ξ)bi + E M [f GP ;2N ]alias (ξ).

(45)

i=0

表 7: CPU time on the computation of the convecting
vortex using SCHGP (NDOF = 303 )f SP ;M (ξ) =

M
∑
i=0

Order of SP M
CPU time

6.

1396[s]

2
246[s]

3

4

198[s]

bi =182[s]

Sec.4. 及び 5. では修正関数に gGa を用い，NC が解の
次数 M ≥ 2N の時に一様流を保持する事を示した．しかし，修正関数に式 (10) の gDG を用いて Sec.5. の数値実験を行うと，表 8 のように M ≥ 2N − 2 の時に一様流を保持する事が分かった．本章では，gDG を修正関数に用いた場合に，gGa より低い解の次数で一様流を保持する
理由を考察する．

−1

Pi (ξ)f GP ;2N (ξ)dξ

(46)

[
]
f C (ξ) = f SP ;M (ξ) + f GP ;2N (−1) − f SP ;M (−1) gL (ξ)[]+ f GP ;2N (+1) − f SP ;M (+1) gR (ξ)(47)(−1)M +1
{PM +1 (ξ) − PM (ξ)}
(48)
2
M +1
(−1)
{PM +1 (−ξ) − PM (−ξ)}
gR (ξ) = gL (−ξ) =2(49)
gL (ξ) =表 8: Freestream preservation using N C metric evaluation with gGa and gDG modiﬁed fuction. The order of ashape fuction is N = 2.

NC (with gGa )∫ 1

E M [f GP ;2N ] は高次モード（2N 次）が低次モード（M次以下）に及ぼすエイリアシング誤差である（補遺 A）．
式 (18) の再構築メトリックを想定し，以下の再構築関数
f C (ξ) を定義する．

NC における修正関数と一様流保持性の関係

Order of SP M

2N + 1
2

NC (with gDG )1 (M = 2N − 3)−2

5.2 × 10

8.8 × 10−3

2 (M = 2N − 2)
3 (M = 2N − 1)

1.9 × 10−2
4.4 × 10−3

4.4 × 10−14
7.6 × 10−14

4 (M = 2N )1.5 × 10−13

1.6 × 10−13

ここで，gL,R は左右の修正 Radau 関数である（Sec. 2.1）．

5
以下の議論は前節で行った再構築メトリックで用いた多次元補間にも拡張出来ると考えられる．

9

Copyright c 2013 by JSFM

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

6.1

M = 2N − 1 の場合
2N −1

GP ;2N

補遺 A-1 を用い，E[fると以下のように変形出来る．
2N−1∑

f C (ξ) =1. 再構築した流束が一様流保持（面積保存則を満足）する為には，SP 及び Modiﬁed part（式 (19)）がそ
れぞれ零となれば十分である事を示した．また，保存量保存性は各セル境界でメトリックの値が一意に定まる事が十分である事を示した．

]alias (ξ) = 0 を代入す

Pi (ξ)bi + P2N (ξ)b2N

i=0
2N
∑

=

Pi (ξ)bi = f GP ;2N (ξ)

(50)

i=0

このように M = 2N − 1 次で再構築した関数 f C (ξ) は，元の 2N 次関数の値 f GP ;2N と一致する．これは再構築
メトリック ξˆxC が，形状関数で定義される厳密なメトリック ξˆxGP ;2N と一致する事に対応する為，SP で面積保存則を満足する．

6.2

M = 2N − 2 の場合

補遺 A-2 を用い，E 2N −2 [f GP ;2N ]alias (ξ) = −(2N −2)P2N −2 (ξ)b2N /2N を代入すると以下のように変形出来る．

f C (ξ) =2N−1∑

Pi (ξ)bi + P2N −2 (ξ)b2N

3. SCHGP を用いると任意の次数で一様流と移流渦が保持される事を示した．一方，一様流は保持するがメトリックの次数が解の次数よりも低い場合（SC を使用した場合），移流渦の解像度は悪くなり，数値解の精度がメトリックの次数で決まる事を示した．

(51)

i=0

4. NC の修正関数に Radau 多項式を使用すると，解の次数が形状関数の 2 倍未満でも一様流を保持する事を数値的に示した．これは，修正関数に Radau 多項式を用いた NC が解析的なメトリック（形状関数の次数の 2 倍の次数を持つ）に一致し，SP 及び Modiﬁed
part が打ち消し合って零となる場合に生じる事を示した．

ここで，以下のルジャンドル多項式の微分漸化式を用いる．ξgn は n 次のガウス点（Pn+1 (ξgn ) = 0）である．
0(2n + 3)Pn+1 (ξ) = Pn+2(ξ) − Pn0 (ξ)∴

0(ξgn ) = Pn0 (ξgn )
Pn+2

(52)
(53)

これにより式 (51) の一階微分の値は 2N − 2 次のガウス点において以下のように書ける．

f 0C (ξg2N −2 ) =

2N
−1
∑

今後は対称保存型メトリックと保存型 FR 法を用いて，移動変形格子上で高次精度の流体計算を行える手法を提案していきたい．

02N −2
Pi0 (ξg2N −2 )bi + P2N)b2N−2 (ξg謝辞
本研究は JSPS 科研費 258793 の助成を受けた．ここに記して謝意を表する．

i=0

=

2N
−1
∑

0
Pi0 (ξg2N −2 )bi + P2N(ξg2N −2 )b2N

補遺 A：補間多項式の次数の違いに伴うエイリアシング

i=0

=

2N
∑

Pi0 (ξg2N −2 )bi = f 0GP ;2N (ξg2N −2 )

i=0

(54)
すなわち，M = 2N − 2 次で再構築した関数 f C (ξ) は，
2N − 2 次のガウス点において一階微分が元の 2N 次関数の一階微分 f 0GP ;2N と一致する．これは面積保存則に現れる再構築メトリック ξˆxC の各方向微分項が，元の ξˆxGP ;2N
等の各方向微分項と一致する事に対応し，SP において一様流を保持する事が分かる．

7.

2. 従来から用いられている NC（非保存型メトリック）は，解の次数が形状関数の次数の 2 倍以上の時にのみ一様流を保持する事を示した．これに対し，任意の次数で一様流を保持する SCHGP（対称保存型高次メトリック）の構築方法を詳細に説明した．なお，
SCHGP はメトリックの次数と解の次数が常に一致するように構築される．一方，解の次数が形状関数の次数よりも高い時に一様流を保持するものの，メ
トリックの次数が形状関数の次数以上にならない SC（対称保存型メトリック）を構築する手法も説明した．

誤差
簡単の為，1 変数関数で説明する．本スキームでは選点多項式を用いて補間を行う為，与えられた N +1 個の離散値から回復される N 次多項式は補間方法によらず一意に定まる．以下では離散値は全て必要な次数のガウス点で定義される（本論では SP に相当）．f N (x) を x ∈ [−1, 1]NN
Nの N 次関数とすると，xN
i = {x0 , x1 , . . . , xN }（N + 1
個の N 次ガウス点）における離散値から構築される N
次の補間多項式 I N [f N ](x) で厳密に表現出来る．

f N (x) = I N [f N ](x) =

結論

N
∑
i=0

3 次元任意形状の静止格子で高次精度保存型 FR 法を用いる場合に一様流を保持する保存型メトリックの評価
方法を詳細に示し，検証計算を行った．本研究で得られた結論は以下である．

NN
φNi (x)f (xi ) =N∑Pi (x)bi ,

i=0

(55)
ここで，φNi (x) は xi における N 次のラグランジュ関数，
Pi (x) と bi は i 次モードのルジャンドル多項式と係数で10

Copyright c 2013 by JSFM

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

式 (61) では，N 次関数の厳密な微分を行い M 次ガウス点で離散値 ∂x f N を得て，M 次で補間多項式を再構築する．一方，式 (62) では M 次ガウス点で離散値 f N を得て，M 次で補間多項式を構築し，M 次関数の厳密な微分を行う．詳細は省略するが，M < N の場合には式 (61)は N − 1 次を M 次で補間して M 次，式 (62) は N 次をM 次で補間し微分して M − 1 次となり最高次に差が生じる為，df a (x) 6= df b (x) である．一方，M ≥ N の場合には df a (x) = df b (x) が成立し，微分と補間は可換である．

ある．bi は N 次のガウス点における f N の離散値と Piから定まる．以下ではルジャンドル多項式を用いて，補間次数の違いに伴う誤差を定量的に見積もる．

A-1: N − 1 次での補間N 個の N − 1 次ガウス点における離散値から N − 1次の多項式 I N −1 [f N ](x) を構築すると，以下のように書ける．
I N −1 [f N ](x) =N−1∑

Pi (x)bi ,

(56)

i=0

N − 1 次ガウス点は N 次のルジャンドル多項式の零点である為，I N −1 [f N ](x) には N 次モードの影響は入らない．
f N (x) と I N −1 [f N ](x) の関係は以下のように書ける．

補遺 C：対称保存型メトリックと面積保存則
以下では N, M は補間多項式の次数を表す．対称保存型メトリックの補間多項式を以下のように定義する（式(33),(35) に相当）．

ξˆxN = ∂ζ I N [TA ] − ∂η I N [TB ]

(63)

(57)

η̂xN = ∂ξ I N [TB ] − ∂ζ I N [TC ](64)

A-2: N − 2 次での補間
N − 1 個の N − 2 次ガウス点における離散値から N − 2次の多項式 I N −2 [f N ](x) を構築すると，以下のように書ける．

ζ̂xN = ∂η I N [TC ] − ∂ξ I N [TA ]

(65)

N

f (x) = I

I

N −2

N

[f ](x) =

N −1

N
−2
∑

N

[f ](x) + PN (x)bN解の次数が M の時，面積保存則 (15) の右辺は以下のよ
うに計算される．

∂ξ I M [ξˆxN ] + ∂η I M [η̂xN ] + ∂ζ I M [ζ̂xN ]Pi (x)bi + E

N −2

N

[f ]alias (x), (58)i=0

E N −2 [f N ]alias (x) = −
N −2

N −1
PN −2 (x)bN ,N(66)= ∂ξ I [∂ζ I [TA ]] − ∂ζ I [∂ξ I [TA ]]
M

N

M

N

+ ∂η I M [∂ξ I N [TB ]] − ∂ξ I M [∂η I N [TB ]](59)+ ∂ζ I M [∂η I N [TC ]] − ∂η I M [∂ζ I N [TC ]]E[f ]alias (x) を，N − 2 次以下のモードで表現出来ない高次モード (N 次) が N − 2 次以下のモードを汚染するエイリアシング誤差と定義する．選点にガウス点を用いるとルジャンドル多項式の性質により，式 (59) のようにエイリアシング誤差を陽に表現出来る．エイリア
シング誤差が汚染するモードの係数は正しい値（十分な次数で補間多項式を作った場合の係数）から外れる．式(59) の場合には N − 2 次のモードの係数が，bN −2 →bN −2 − (N − 1)/N と汚染される事が分かる．N − 2 次の補間多項式 I N −2 [f N ](x) を用いると，f N (x) は以下のように回復される．

(67)

N

式 (67) の各行が 0 となれば，面積保存則は満足される．
代表して下線部を書き下すと以下のようになる．

∂ξ I M [∂ζ I N [TA ]] − ∂ζ I M [∂ξ I N [TA ]]MM
dφM
l (ξ)φm (η)φn (ζ)

l,m,n=0N∑NN
φNi (ξl )φj (ηm )dφk (ζn )TA;i,j,k

i,j,k=0−

f N (x) = I N −2 [f N ](x) − E N −2 [f N ]alias (x)+ PN −1 (x)bN −1 + PN (x)bN

M
∑

=

(68)

M
∑

M
M
φM
l (ξ)φm (η)dφn (ζ)

l,m,n=0
N
∑

(60)

N
N
dφN
i (ξl )φj (ηm )φk (ζn )TA;i,j,k(69)

i,j,k=0

なお，任意の次数の補間多項式を用いた際のエイリアシ
ング誤差の議論は本稿の主旨から外れるので省略するが，今後別の論文で紹介したい．
補遺 B：微分を含む補間多項式の次数の違いに伴う非可
換性

(61)

df b (x) = ∂x I M [f N ](x).

(62)

M ∑
N
∑

N
φM
m (η)φj (ηm )m=0 j=0[ M N∑∑
n=0 k=0|先と同様に，N 次関数 f N (x) を考える．以下の 2 つの補間多項式を考える．

df a (x) = I M [∂x f N ](x),=NφMn (ζ)dφk (ζn ){z}(ii-a){ M N∑∑
l=0 i=0

|

}]

N
dφM
l (ξ)φi (ξl ) TA;i,j,k

{z

}

(i-a)

11

Copyright c 2013 by JSFM

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

−

M ∑
N
∑

(8) Deng, X., Mao, M., Tu, G., Liu, H., and Zhang, H.,“Geometric Conservation Law and Applications to
High-Order Finite Diﬀerence Schemes with Stationary Grids,” Journal of Computational Physics,
Vol. 230, No. 4, 2011, pp. 1100–1115.

N
φMm (η)φj (ηm )m=0 j=0[ M N∑∑
n=0 k=0|N
dφMn (ζ)φk (ζn ){z}(ii-b){ M N∑∑
l=0 i=0|(9) Deng, X., Min, Y., Mao, M., Liu, H., Tu, G., andZhang, H., “Further studies on Geometric Conservation Law and applications to high-order ﬁnite difference schemes with stationary grids,” Journal ofComputational Physics, Vol. 239, 2013, pp. 90–111.

}]

N
φM
l (ξ)dφi (ξl ) TA;i,j,k

{z

(70)

}

(i-b)

式 (70) が 0 となる為には，(i-a) = (i-b) かつ (ii-a) = (ii-b)が十分条件である．補遺 B で説明したように M ≥ N であれば微分と補間が可換となり，(i-a) = (i-b) かつ (ii-a) =(ii-b) が成立する．逆に，M < N ならば (i-a) 6= (i-b) かつ (ii-a) 6= (ii-b) である為，式 (70) は 0 とならない．すなわち，式 (63) の形で書かれる対称保存型メトリックが面積保存則を満足するのは M ≥ N の時である．

(10) Abe, Y., Nonomura, T., Iizuka, N., and Fujii, K.,“Geometric interpretations and spatial symmetry
property of metrics in the conservative form forhigh-order ﬁnite-diﬀerence schemes on moving anddeforming grids,” Submitted to Journal of Computational Physics, 2013.
(11) Nonomura, T., Iizuka, N., and Fujii, K.,“Freestream and vortex preservation propertiesof high-order WENO and WCNS on curvilinear
grids,” Computers & Fluids, Vol. 39, No. 2, 2010,
pp. 197–214.

参考文献(1) Vinokur, M. and Yee, H., “Extension of Eﬃcient Low Dissipation High Order Schemes for 3-D
Curvilinear Moving Grids,” Frontiers of Computational Fluid Dynamics, 2002, pp. 129–164.

(12) Nonomura, T., Terakado, D., Abe, Y., and Fujii, K., “A Technique for Freestream Preservationof Finite-Diﬀerence WENO on Curvilinear Grid,”
Submitted to Journal of Computational Physics,
2013.

(2) Abe, Y., Iizuka, N., Nonomura, T., and Fujii, K.,“Conservative metric evaluation for high-order ﬁnite diﬀerence schemes with the GCL identities onmoving and deforming grids,” Journal of Computational Physics, Vol. 232, 2013, pp. 14–21.

(13) Kopriva, D., “Metric Identities and the Discontinuous Spectral Element Method on Curvilinear
Meshes,” Journal of Scientific Computing, Vol. 26,
2006, pp. 301–326.

(3) Wang, Z. J., “High-Order CFD Methods: Current
Status and Perspective,” International Journal for
Numerical Methods in Fluids, Vol. 00, 2012, pp. 1–
42.

(14) Kim, J. W. and Lee, D. J., “Characteristic Interface Conditions for Multiblock High-Order Computation on Singular Structured Grid,” AIAA Journal , Vol. 41, 2003, pp. 2341–2348.

(4) Huynh, H. T., “A Flux Reconstruction Approach
to High-Order Schemes Including Discontinuous
Galerkin Methods,” AIAA 2007-4079 , 2007.

(15) Sumi, T., Kurotaki, T., and Hiyama, J., “Generalized Characteristic Interface Conditions forHigh-order Multi-block Computation,” International Journal of Computational Fluid Dynamics,
Vol. 21, No. 1, 2007, pp. 9–10.

(5) Haga, T. and Kawai, S., “Toward Accurate Simulation of Shockwave-Turbulence Interaction on Unstructured Meshes: A Coupling of High- Order FR
and LAD Schemes,” AIAA 2013-3065 , 2013.
(6) Thomas, P. D. and Lombard, C. K., “Geometric Conservation Law and its Application to Flow
Computations on Moving Grids,” AIAA Journal ,Vol. 17, No. 10, 1979, pp. 1030–1037.

(16) Kawai, S. and Lele, S., “Localized artiﬁcial diﬀusivity scheme for discontinuity capturing on curvilinear meshes,” Journal of Computational Physics,
Vol. 227, 2008, pp. 9498–9526.

(7) Visbal, M. R. and Gaitonde, D. V., “On the Use ofHigher-Order Finite-Diﬀerence Schemes on Curvilinear and Deforming Meshes,” Journal of Computational Physics, Vol. 181, No. 1, 2002, pp. 155–
185.
12

Copyright c 2013 by JSFM