第 28 回数値流体力学シンポジウム
C06-4

３次元移動変形格子での高次精度保存型 FR 法における一様流保持
Freestream preservation of a high-order conservative FR scheme on moving and deforming grids
阿部圭晃，東大院, 相模原市中央区由野台 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 flux reconstruction (FR) scheme on a three-dimensionally-curved moving and deforming
grid. In our consistent symmetric-conservative metric evaluation, the metrics is appropriately computed to avoid
the freestream preservation errors for any couple of the order of a solution polynomial and high-order shapefunctions. We computed a convecting vortex using compressible Euler equations on three-dimensionally deforming
wavy meshes, where the freestream preservation and global conservation properties were demonstrated.

1.

はじめに

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

いても，セルの形状を定義する形状関数が高次である場
合には同様の工夫が必要であるが(10)(11) ，移動変形を伴う格子への適用例はこれまでに存在しない．本研究では，
移動変形を伴う格子における保存型メトリックを FR 法に拡張し，保存量保存性と一様流保持性を両立する手法を提案する．

2.

保存型 Flux Reconstruction 法

2.1

支配方程式の離散化

計算領域は各時刻において空間六面体セルで分割されているものとする．デカルト座標系 (x, y, z, t) で表される物理空間の全てのセルを，それぞれ一般座標 (ξ, η, ζ, τ )で表される計算空間の標準セルに変換する．各座標系の間には以下の関係を仮定する．

ξ = ξ(x, y, z, t), η = η(x, y, z, t), ζ = ζ(x, y, z, t),(1)x = x(ξ, η, ζ, τ ), y = y(ξ, η, ζ, τ ), z = z(ξ, η, ζ, τ ), (2)t = τ,(3)

物理空間の第 n セルの形状を表す為に，一般座標の空間
各方向毎に N + 1 個ずつ，時間方向に T + 1 個の Grid
point (以下 GP) を定義し，それらの座標 rn;i,j,k,l に対応
する空間方向 N 次，時間方向 T 次の形状関数を用いる（図 1）．なお，説明の都合上時間方向の形状関数を導入するが，本稿では差分法に基づく通常の時間積分を用いる1 ．また本稿では，簡便化の為に以下のような補間多項式の表現を用いる．

ψ GP ;m,n (ξ, η, ζ, τ ) := I GP ;m,n [ψ(ξ, η, ζ, τ )],(4)

ψ SP ;m,n (ξ, η, ζ, τ ) := I SP ;m,n [ψ(ξ, η, ζ, τ )].

(5)1
後述するが，時間方向の解の定義点は GP に一致する．時間方向には，時間積分に用いる全ての離散時間ステップ（例えば，Runge-Kutta 法では全てのsubstep）を含むようにセルを仮定する．

Copyright c 2014 by JSFM

第 28 回数値流体力学シンポジウム
C06-4

式 (4),(5) はそれぞれ，任意の多項式関数 ψ を GP・SP における離散値から再構築した空間 m 次・時間 n 次の補間
多項式 ψ GP ;m,n , ψ SP ;m,n である．この表現を用いると，
セル形状は空間 N 次，時間 T 次の補間多項式を用いて以下のように書ける．





ρV



ρW







 ρuV + ηx p  ρuW + ζx p 







F̂ = J  ρvV + ηy p  , Ĝ = J  ρvW + ζy p  , ρwV + ηz p  ρwW + ζz p (e + p)V − η̂t p(e + p)W − ζ̂t p

rn (ξ, η, ζ, τ ) = I GP ;N,T [rn ]
=

N
T
∑
∑

U = ξˆt + ξx u + ξy v + ξz w,

N
N
N
T
Mn;i
(ξ)Mn;j(η)Mn;k(ζ)Mn;l(τ )rn;i,j,k,lV = η̂t + ηx u + ηy v + ηz w,

i,j,k=0 l=0(6)

W = ζ̂t + ζx u + ζy v + ζz w.
式 (7) を解く為に，各セルに内点（内部自由度）を導入
し解を多項式近似する．一般座標の空間各方向毎に K + 1
個，時間方向には T + 1 個（GP の時間方向離散点と一致）の Solution point （以下 SP：図 1 の緑色の点）を定義し，各方向のラグランジュ多項式のテンソル積として保存量 Q̂ の第 n セル内の分布を与える．

τl

;K,T
Q̂SP
:= I SP ;K,T [Q̂n ]
n

=

τl+1

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

座標変換のヤコビアン J := |∂(x, y, z)/∂(ξ, η, ζ)| とメ
トリック ξˆx = ξx J, ξˆt = ξt J 等はこれらの微分を元に計算出来る．なお，図 1 は N = 1, T = 1 の場合の 1 セルを示しており，GP は青色のセル頂点にのみ存在する．N ≥ 2の場合にはセル境界上にも GP が定義されるが，実際に計算格子を準備する際には各時刻で頂点の GP のみを定
義すればよい（空間的なセル境界上の GP は，計算空間
で等間隔に置けば良い）．なお今回の論文では格子の移動変形を解析式で定義するが，一般の形状周りに格子を作る時には高次の GP も何らかの方法で定義する必要がある．計算空間における圧縮性 Euler 方程式は次のように書ける2 ．

2

T
K
K
φK
n;p (ξ)φn;q (η)φn;r (ζ)φn;l (τ )Q̂n;p,q,r,l (9)p,q,r=0 l=0

図 1: Deﬁnition of grid, solution and ﬂux points in eachcell ( spatially 2D case is shown). The order of a shapefunction is N = 1 in space and T = 1 in time; the orderof a solution polinomial is K = 2 in space and T = 1 intime, respectively.

∂τ Q̂ + ∂ξ Ê + ∂η F̂ + ∂ζ Ĝ = 0,
 


ρ
ρU
 


 ρu 
 ρuU + ξx p 
 





Q̂ = J  ρv  , Ê = J  ρvU + ξy p  ,
 


ρw
 ρwU + ξz p e(e + p)U − ξˆt p

M
T
∑
∑

ÊnSP ;K (ξ) := I SP ;K [ÊnSP ](ξ)
=

K
∑

SP
φK
n;p (ξ)Ên;p,q,r,l ,

(10)

p=0

∂ξ ÊnSP ;K (ξ) := ∂ξ I SP ;K [ÊnSP ](ξ)
=

K
∑

SP
dφK
n;p (ξ)Ên;p,q,r,l ,

(11)

p=0

φK
n;p は以下で決まる K 次ラグランジュ多項式であり，
K
dφKn;p = dφn;p /dξ とする．

(7)

φKn;p (ξ) =K∏
m=0,m6=p(8)

ξ − ξm.
ξp − ξm(12)

この流束 ÊnSP ;K は隣接セルの情報を反映せず，セル境界で不連続な分布となる．そこで，セル境界で共通の値を採るような分布 ÊnC を再構築する．まず，セル境界に各方向の SP を通る直線とセルを囲む面の交点に Flux point（以下 FP：図 1 の赤色の点）を置く．式 (9) を用いて境界

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

2

Copyright c 2014 by JSFM

第 28 回数値流体力学シンポジウム
C06-4の両側のセルから FP における解 Q̂L , Q̂R を外挿し，近似リーマン解法を用いて Ê com を風上的に求める．そして，FP で Ê com を通るように元の K 次の ÊnSP ;K を修正
し ÊnC を得る.
][
com
− ÊnSP (−1) gL (ξ)
ÊnC (ξ) = ÊnSP (ξ) + Ên−1/2[]
com+ Ên+1/2− ÊnSP (+1) gR (ξ) (13)
ここで，左右のセル境界は n ∓ 1/2 で示す．gL (ξ) は左端（ξ = −1）で 1，右端（ξ = +1）で 0 となる K + 1 次多項式であり，修正関数(2) と呼ばれる．本稿では以下で定義される関数を採用し，gL = gGa,K+1 とする．

K
K −1
RR,K +
RR,K−1 ,2K − 1
2K − 1(−1)K
RR,K = gDG,K =(PK − PK−1 ),2

gGa,K =

(14)
(15)

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

なお，FP において共通のメトリックを使用せずに一様流を保持するスキームを構築する事は可能であるが，その場合一般には保存量保存性が満たされない事を強調しておきたい．各セル内で独立に（局所的に）議論される
一様流保持と，隣り合うセルの持つ値まで考慮して議論される保存量保存性の概念を両立する為には特別な工夫が必要である．

FR 法における一様流保持の条件

3.

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

3.1

式 (13) で得られる流束 ÊnC は境界を含むセル内で連続な
M 次の多項式であり，式 (16) は満足する．SP としてガ
ウス点を用いている為，式 (13) の ÊnSP (±1) は式 (11) のÊnSP により外挿された値である．この時，保存量 Q を外挿して境界の流束を計算すると内部の近似多項式との不整合を生じ，局所保存性が失われる事に注意する．
次に，計算領域内での保存性（保存量保存性）を議論する．任意の閉領域 ∀ Ω における保存量保存性が満たされるには次式を満たす必要がある．
∑ ∫∂ξ ÊnC dξ = Ê∂com− Ê∂com
.
∀ Ω;n
∀ Ω;n
max +1/2
min −1/2n∈∀ Ωcom
もし各セル境界（FP）の Ên±1/2について，隣り合うセル同士で同じ値を共有すれば式 (18) 左辺の総和は閉領域∀Ω 境界の値のみが残る．ゆえに，各セル境界（FP）でcom
Ên±1/2の値が一意に定まる事は保存量保存性が成り立つ
為の十分条件である．これは直交等間隔格子では自明に
com成り立つ．しかし格子形状が曲線の場合には，Ên±1/2に含まれるメトリックもセル同士で共有する FP で同じ値を持たなければならず自明でない．本研究では，FP でのメトリックを構築する際に隣り合うセル同士で同じ離散値から補間して与える(10) ．例えば，ξx J = yη zζ − yζ zη
を n + 1/2 番目の FP で定義する際には，n, n + 1 番目の両セルで共有する GP の値（y, z ）のみを用いて補間関数を構築すれば良い．このように隣り合うセルで GP を共有させる為に，各方向に分布させる GP の端点は必ずセルの境界に位置していなければならない．本稿のメト
リックはこれらの性質を満たすように構築しており，保存量保存性が保証されている．

一様流保持と体積保存則

一様流計算時の保存量ベクトル Q̂(ξ)，再構築流束Ê (ξ) を書き下す．なお以下は全て時刻 τ = τl における第 n セル内の議論であり，下付き添字の n, l は一部省略する．
CQ̂(ξ) = Quni J SP (ξ)
(19)
[
]
com
Ê C (ξ) = Ê SP (ξ) + Ên−1/2− Ê SP (−1) gL (ξ)[]
com+ Ên+1/2− Ê SP (+1) gR (ξ)Ωn= Quni ξˆtC (ξ) + E uni ξˆxC (ξ)(17)+ F uni ξˆyC (ξ) + Guni ξˆzC (ξ)∀

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

ここで，Quni , E uni , F uni , Guni は一様流の値，ξˆxC ={ξx J}C 等 は 後 述 す る 再 構 築 メ ト リック で あ る ．式(19),(20) を支配方程式 (7) に代入すると，以下の体積保存則（volume conservation law: VCL）と面積保存則（sur-n∈∀ Ω. (18)− Ê∂com= Ê∂com
∀ Ω;n
∀ Ω;n
min −1/2
max +1/2
3

Copyright c 2014 by JSFM

第 28 回数値流体力学シンポジウム
C06-4における離散値を計算する．

face closure law: SCL）を得る．
∂τ J SP + ∂ξ ξˆtC + ∂η η̂tC + ∂ζ ζ̂tC = 0,(21)

J BS;3N −1,3T∂ξ ξˆxC + ∂η η̂xC + ∂ζ ζ̂xC = 0,(22)

= (∂ξ I GP ;N,T [r] × ∂η I GP ;N,T [r]) · ∂ζ I GP ;N,T [r], (29)∂ξ ξˆyC + ∂η η̂yC + ∂ζ ζ̂yC = 0,(23)

ζ̂tBS;3N,3T −1∂ξ ξˆzC + ∂η η̂zC + ∂ζ ζ̂zC = 0,(24)

= −(∂ζ I GP ;N,T [r] × ∂η I GP ;N,T [r]) · ∂τ I GP ;N,T [r],(30)

これらを併せて幾何学的保存則（geometric conservation
law: GCL）と呼び，これら全てを満たす事が一様流保持の十分条件となる．以下では体積保存則（VCL）(21) に注目する．面積保存則（SCL）(22)–(24) を満たす空間メ
トリックには過去研究で提案したもの(10) を用いる．

上付き添字が示すように，J BS は空間 3N − 1 次，時間
3T 次の補間多項式として定義される．同様に，ζ̂tBS は空間最大 3N 次，時間 3T − 1 次の補間多項式となる．これ
らはセル形状を厳密に表現するヤコビアンとメトリック
であり，充分な次数（空間 3N 次，時間 3T 次）の補間多項式に基づき微分を行えば，連鎖則（Chain rule）と多重微分の可換性（commutativity）が満たされ(10) ，VCL を満たす．しかし，実際には SP，FP における J BS と ζ̂tBSの離散値を用いて保存量ベクトルと流束を計算し，流束の空間微分及び時間積分は SP の次数（空間 K ，時間 T ）の多項式に基づいて行う．ゆえに，SP part（VCL の一部）の各項に現れる微分項は J BS と ζ̂tBS を解析的に微分したものではなく，J BS と ζ̂tBS を SP の次数に落とした多項式を解析的に微分したものである．このように，J BSと ζ̂tBS を SP の次数に落とした多項式が，先に導入した
J SP と ζ̂tSP に相当し，以下のように計算される：

体積保存則の成立条件
先に示したメトリック ξˆtC 等を以下では再構築メトリッ
クと呼ぶ．式 (20) に従うと，ξˆC は

3.2

t

[
]
ξˆtC (ξ) = ξˆtSP (ξ) + ξˆtF P (−1) − ξˆtSP (−1) gL (ξ)][+ ξˆtF P (+1) − ξˆtSP (+1) gR (ξ) (25)と書ける．ここで，ξˆtSP は SP で，ξˆxF P は FP で用いられるメトリックである．引数の ±1 は，各セルの左右のFP を示す．式 (25) より体積保存則 (21) を書き下すと以下の通りである．

J SP ;K,T = I SP ;K,T [J BS;3N −1,3T ]式 (21) 左辺 ⇐⇒∂τ J SP + ∂ξ ξˆtSP + ∂η η̂tSP + ∂ζ ζ̂tSP
SP part[]FP
SP+ ξˆt (−1) − ξˆt (−1) dgL (ξ)/dξ + . . .

=

K
T
∑
∑

BS;3N −1,3T
T
K
K
,
φK
p (ξ)φq (η)φr (ζ)φs (τ )Jp,q,r,sp,q,r=0 s=0(31)

correction partζ̂tSP ;K,T = I SP ;K,T [ζ̂tBS;3N,3T −1 ](26)
式 (26) の下線部（SP 及び correction part）が SP において 0 となれば体積保存則は満たされ，一様流保持が成
立する．

4.

=

BS;3N,3T −1
K
K
T
φK
,
p (ξ)φq (η)φr (ζ)φs (τ )ζ̂t;p,q,r,sp,q,r=0 s=0(32)

時間メトリックとヤコビアンの構築方法
上に示す J SP ;K,T , ζ̂tSP ;K,T では空間方向に 3N → K ，時間方向に 3T → T の近似を行っている．ここでは，空間方向にはセル形状を定義する多項式の次数と解を表現する多項式の次数を独立に定める（N 6= K を許す）一方
で，時間方向にそのような特殊な取り扱いをしない事に注意する．この為，K ≥ 3N と選べば，空間方向に厳密な表現となる（3N → K = 3N ）一方3 で，時間方向には必ず低次近似（3T → T ）を伴う．この為，NC を用いた
場合には J SP ;K,T , ζ̂tSP ;K,T は時間方向に必ず打ち切り誤差やエイリアシング誤差を含み，微分の連鎖則と可換性が成り立たず式 (26) の SP part は 0 とならない．ゆえに，
NC を用いると一般には移動変形格子上で一様流を保持しない．なお，先に述べたように時間方向にもセル形状の定義多項式の次数を解の多項式の次数と独立に定める

以下では時刻 τ = τs における第 n セル内の SP・FP におけるメトリックとヤコビアンの構築方法を述べる．以下では，GP を i, j, k, l，SP を p, q, r, s のインデックスで表す（先にも述べたが，本稿では時間方向の SP と GP は一致するため，l と s の示す離散点は同一である）．

4.1

K
T
∑
∑

Nonconservative metric (NC)

一般に，時間メトリックとヤコビアンは以下のように定義される事が多い．これらの表式を非保存型（nonconservative metric: NC）と呼ぶ．

(rξ × rη ) · rζ ,

ξˆt = −(rη × rζ ) · rτ ,(27)

η̂t = −(rζ × rξ ) · rτ ,

ζ̂t = −(rξ × rη ) · rτ .

(28)

J=

以下では，J と ξˆt に着目する．まず GP の値を用い，以下のように J BS , ζ̂tBS なる補間多項式を構築し，SP，FP

K ≥ 2N の時に NC が SCL を満足する例に相当する（Eq.(44) in Sec.5.2
）

3
(10)

4

Copyright c 2014 by JSFM

第 28 回数値流体力学シンポジウム
C06-4

ような定式化を行えば，この限りではない．また本稿のNC の定式化を用いても，セル形状が時間的に変化しな
い移動格子を用いる（格子が剛体運動する：ヤコビアン
が時間的に変化しない）場合には，空間方向に十分な次数（K ≥ 3N ）を保てば一様流を保持する（補遺参照）．

4.2 Symmetric conservative metrics
一様流保持を満たす為に，以下のような対称保存型（symmetric conservative metrics）の時間メトリックとヤコビアンを導入する．これは，NC（式 (32)）を解析的に書き換えたものである．
J = [{(rξ × r)η · r + (r × rη )ξ · r}ζ

(33)+{(rζ × r)ξ · r + (r × rξ )ζ · r}η

(34)+{(rη × r)ζ · r + (r × rζ )η · r}ξ ]/6(35)ξˆt = [−{(rη × r)ζ · r + (r × rζ )η · r}τ

(36)−{(rτ × r)η · r + (r × rη )τ · r}ζ

(37)−{(rζ × r)τ · r + (r × rτ )ζ · r}η ]/6(38)η̂t = [{(rξ × r)ζ · r + (r × rζ )ξ · r}τ

(39)+{(rζ × r)τ · r + (r × rτ )ζ · r}ξ

(40)+{(rτ × r)ξ · r + (r × rξ )τ · r}ζ ]/6(41)

ζ̂t = [−{(rξ × r)η · r + (r × rη )ξ · r}τ

(42)−{(rτ × r)ξ · r + (r × rξ )τ · r}η

(43)−{(rη × r)τ · r + (r × rτ )η · r}ξ ]/6(44)

以上より，FR 法で移動変形格子を用いる場合に一様流保持を満たす為の要点は，(1) 対称保存型に書かれた時間メトリックとヤコビアンの被微分項（式 (47)，(54) の G ）を同一の多項式で表す(2) 支配方程式の時間・空間微分の近似次数（T, K ）を，構築した時間メトリックやヤコビアンの補間多項式の次数以上にする
事の 2 点である．(1),(2) に基づき，本稿では standard
symmetric conservative metrics (SSC) と呼ばれる時間
メトリックとヤコビアンの定式化を提案する．本稿ではこれに加えて，(3) 精度を保証する為に，時間メトリックとヤコビアンの補間多項式の次数を解の近似次数（T, K ）に一致させる4という要請を導入する（2 点目の要請をより厳しいものに置き換えたと考えて良い）．(1),(3) に基づき，consistentsymmetric conservative metrics (CSC) をと呼ばれる定式化も提案する．

4.2.1 SSC ヤコビアン J BS と時間メトリック ζ̂tBS は
GP の値を用いて以下のように計算される．簡単のため，行 (33) と行 (42) のみ示す．
J BS;N,T = [∂ζ I GP ;N,T [G GP ;2N,2T ] + · · · ,(45)ζ̂tBS;N,T = [∂τ I GP ;N,T [G GP ;2N,2T ] + · · · ,(46)

G GP ;2N,2T は GP での離散値を用いた空間最大 2N 次，時間 2T 次の補間多項式であり，以下のように定義される：

さて，NC の時間メトリックとヤコビアンが VCL を満たすには微分の連鎖則と可換性の両方が必要であり，メトリッ
クとヤコビアンの補間多項式を十分な次数（空間 3N − 1
次，時間 3T 次）で表現する必要があった．しかし，対称保存型の場合は多重微分の可換性のみを必要とする為，比較的容易に VCL を満たす事が出来る．例えば J の行 (33)を τ で微分した項：{(rξ ×r)η ·r +(r ×rη )ξ ·r}ζτ と，ζ̂t の行 (42) を ζ で微分した項：{(rξ × r)η · r + (r × rη )ξ · r}τ ζは解析的に打ち消し合う．被微分項に相当する補間多項式は一致していなければならないが，この時，微分の連鎖則は用いておらず，J の最も外側に現れる空間微分（ζ ）と支配方程式の時間微分（τ ）の可換性，ζ̂t の最も外側に現れる時間微分（τ ）と支配方程式（流束）の空間微分（ζ ）の可換性のみに基づく．FR 法では支配方程式の時
間・空間微分の離散近似の次数が，ζ̂t や J の最も外側の微分を表現する多項式の次数（ζ̂t と J の補間多項式の次数に相当）以上であれば，この可換性が成立し（異なる次数の補間多項式を用いた微分の可換性については参考
文献 (10) の Appendix. C を参照），対称保存型を用いてVCL を満足し一様流を保持する事が出来る．

5;2N,2T
G GP ;2N,2T = ∂η I GP ;N,T [T GP] · I GP ;N,T [r]A;2N,2T+ ∂ξ I GP ;N,T [T GP] · I GP ;N,T [r], (47)B
;2N,2T
T GP= ∂ξ I GP ;N,T [r] × I GP ;N,T [r],A

(48)

;2N,2T= I GP ;N,T [r] × ∂η I GP ;N,T [r]
T GP
B

(49)

このように構築した J BS;N,T , ζ̂tBS;N,T を用いて SP,FP における離散値を計算する．NC の節で説明したように，実際の支配方程式に現れる空間・時間微分は SP の次数に基づく為，J SP , ζ̂tSP は以下のように計算される．

J SP ;K,T = I SP ;K,T [J BS;N,T ],(50)

ζ̂tSP ;K,T = I SP ;K,T [ζ̂tBS;N,T ](51)

対称保存型メトリックを FR で用いる場合に重要な点は，
G の補間多項式を共通に持つ事である（要点 (1)）．これにより，G の離散値を基に再構築・微分される J BS とζ̂tBS の補間多項式が打ち切り/エイリアシング誤差を含むとしても，支配方程式の空間微分の次数 K が J BS と4

空間メトリックに関する結果を基にしている(10)Copyright c 2014 by JSFM

第 28 回数値流体力学シンポジウム
C06-4

ζ̂tBS の補間多項式の次数以上であれば微分の可換性に基づき VCL の SP part は 0 となる（要点 (2)）．ゆえに本
節で導入した SSC では，K ≥ N とすれば VCL の SP
part を満たし（この時，ζ̂tSP ;K,T と ζ̂tBS;N,T は一致する
為，correction part も 0 となる），一様流保持が成り立つ．
しかし，SSC は，式 (47) において G GP ;2N,2T を計算する際，空間の次数を 2N → N とする低次近似を行う．
また，式 (45),(46) において再び空間の次数を 2N → Nと低次近似する．ゆえに，解の空間次数 K を高くしてもJ SP ，ζ̂tSP の次数は N 次より大きくならない．この事は，数値解の空間精度にも影響し，解の空間次数を K > N
としても数値解は N 次に制限される事が予想される(10) ．
すなわち，SSC は要点 (3) を満たさない．
4.2.2 CSC CSC では，ヤコビアンと時間メトリックの空間次数を，SP の空間次数 K に一致させる事を目指す．
その為，空間 K 次の補間多項式を構築する為，空間各方向
毎に K +1 個，時間方向は T +1 個の consistent grid point（以下 CGP）を定義する（時間方向には SP,GP,CGP 共に同じ座標を取る）．CGP からは SP と同じ次数の補間多項式を構築出来るが，CGP と SP が一致する必要はな
い．ヤコビアン J BS と時間メトリック ζ̂tBS は CGP の値を用いて以下のように計算される．簡単のため，行 (33)と行 (42) のみ示す．
J BS;K,T = [∂ζ I CGP ;K,T [G CGP ;K+N,2T ] + · · · ,(52)ζ̂tBS;K,T = [∂τ I CGP ;K,T [G CGP ;K+N,2T ] + · · · ,(53)

G CGP ;K+N,2T は CGP での離散値を用いた空間最大 K +N 次，時間 2T 次の補間多項式であり，以下のように定義される：

特に，K ≥ 3N の場合には空間方向には厳密な時間メト
リックとヤコビアンを回復する5 ．

4.3

NC,SSC,CSC のまとめ

表 1 に，移動変形格子において時間メトリックとヤコビアンに NC,SSC,CSC を用いた場合の一様流保持性をまとめる．

表 1: Expected characteristics of a freestream preservation using each metric evaluation (◦ satisﬁes afreestream preservation; • does not satisfy a freestream
preservation). K is the order of the solution polynomial;
N is the order of the shape function.

Order of SP K

NC

SSC

CSC

K<N
N ≤K

•
•

•
◦

◦
◦

検証計算

5.

前節で定義した NC,SSC,CSC を用いて，一様流保持と移流渦の検証計算を行う．計算対象は圧縮性非粘性流れとし，支配方程式には音速と一様流密度で無次元化された式 (7) を用いる．気体は理想気体とし，比熱比は γ = 1.4とする．時間積分は 2 段 2 次 Runge-Kutta 法を用いる(12)．本計算には，−5 ≤ x, y, z ≤ 5 なる立方体内の直交格子 x0 , y0 , z0 を式 (57)-(59) で 3 次元的に歪めた計算格子を用意し，図 2 に示す．

xi,j,k (τ ) = x0;i,j,k][nπ(l − 1)∆z0nπ(k − 1)∆y0sin, (57)+ ∆x0 α sin
L0L0;2N,2T] · I GP ;N,T [r]G CGP ;K+N,2T = ∂η I CGP ;K,T [T GPAyi,j,k (τ ) = y0;i,j,k[]nπ(j − 1)∆x0nπ(l − 1)∆z0+ ∆y0 α sin
sin
,
L0
L0

;2N,2T] · I GP ;N,T [r], (54)+ ∂ξ I CGP ;K,T [T GP
B

このように構築した J BS;K,T , ζ̂tBS;K,T を用いて SP,FP における離散値を計算する．実際の支配方程式に現れる空間・時間微分は SP の次数に基づく為，J SP , ζ̂tSP は以下のように計算される．

zi,j,k (τ ) = z0;i,j,k[]nπ(k − 1)∆y0nπ(j − 1)∆x0
sin, (59)+ ∆z0 α sin
L0L0
α = A sin(2πωτ ).

J SP ;K,T = I SP ;K,T [J BS;K,T ],(55)

ζ̂tSP ;K,T = I SP ;K,T [ζ̂tBS;K,T ]

(56)

(58)

SSC と同様に，G の補間多項式は J と ζ̂t で共通なので（要点 (1)），微分の可換性に基づいて VCL の SP part は0 となる（要点 (2)）．注意すべきは，式 (55),(56) において低次近似が行われない為，任意の K, N の組み合わせで一様流保持が満たされる事である．これは，J BS と ζ̂tBSの補間多項式の空間次数を SP と同じ K 次で構築した為である．その結果，解の空間次数 K を高くすれば時間メ
トリックとヤコビアンの補間多項式の空間次数も高くなり，解の空間精度が保たれると考えられる（要点 (3)）．
6(60)

n = 4, A = 0.4, L0 = 10, ω = 0.3 とし，∆x0 = ∆y0 =∆z0 はセル頂点間の格子幅とする．流れ場には周期性を仮定し，計算領域の境界では各方向に周期境界条件を課した．
5.1

一様流保持性

初期解 u0 = (u0 , v0 , w0 ), p0 , ρ0 として，x 方向の一様流を以下のように与える．

u0 = 0.10,

v0 = w0 = 0.0,

p0 /γ = ρ0 = 1.0,

(61)

5
T A の次数が 2N 次である為，補間多項式 (54) で定義される G の空間次
数は，（K によらず）3N 次以下である．

Copyright c 2014 by JSFM

第 28 回数値流体力学シンポジウム
C06-4

たすように構築している．実際，表 4 でも保存量保存性の誤差が丸め誤差に抑えられている．

表 4: Global conservation error (hei) evaluated in each
metric. The order of the shape function is N = 2.
Order of SP : K
1
2
34

図 2: Computational grid (the number of cells is
5 × 5 × 5 = 125, and GP = 2 is adopted: total number
of GP including high-order inner grid point is 3375.)
格子点数は各方向に Nmax = 10（自由度は解の次数によって異なる）とした．形状関数の精度が N = 2 の時，10 ス
テップ計算時の流れ場の L∞ ノルム（v, w の最大誤差）を表 2 に示す．表 2 から，NC は SP の次数によらず一

表 2:
L∞ norm of a freestream error(max{|v − v0 |/u0 , |w − w0 |/u0 }). The order of ashape fuction is N = 2.
Order of SP: K
1
2
3
4

NC
4.6 × 10−3
3.1 × 10−3
2.1 × 10−3
1.1 × 10−8SSC
2.0 × 10−3
2.5 × 10−16
9.3 × 10−16
2.2 × 10−14

表 3:
L∞ norm of a freestream error(max{|v − v0 |/u0 , |w − w0 |/u0 }). The order of ashape fuction is N = 1.
Order of SP: K
1
2
3
4

NC
1.9 × 10−3
2.8 × 10−9
9.3 × 10−9
1.8 × 10−8SSC
8.3 × 10−17
9.7 × 10−16
1.5 × 10−16
3.5 × 10−16

CSC
8.3 × 10−17
2.8 × 10−16
4.2 × 10−15
1.9 × 10−14

5.2 保存量保存性
先の一様流に対し，保存量保存性の誤差を表 4 に示す．
物理量 φ に対する保存量保存性の誤差は hφi で定義される（VR は計算領域全体である）
：
/∫
∫
hφi '
(φ − φ0 )dV
φ0 dV ,
(62)VR

u = U∞ − C

VR

領域の積分は Gauss 積分を用いる．Sec.2.2 で述べたように，本稿で提案するメトリックは全て保存量保存性を満

7CSC
3.0 × 10−15
9.8 × 10−15
2.4 × 10−15
2.6 × 10−14(y − yc )
−r2
ρC 2
exp
,
p
=
p
−
exp(−r2 ),
∞
Rc2
2
2Rc2
(63)

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

C = U∞ ,様流保持しない事が分かる．また，SSC は K ≥ N の時に，CSC は任意の K, N の組み合わせで一様流を保持する．これらは事前に示した表 1 と一致する．なお，N = 1の時には表 3 の通りであり，この場合の一様流保持性も表 1 と一致する．

SSC
3.7 × 10−15
9.8 × 10−15
5.6 × 10−15
2.6 × 10−14

5.3 移流渦
次に，2 次元渦(13) を移流させた計算結果を以下に示す．
移流渦は以下のように定義した．

v=C

CSC
6.9 × 10−17
2.9 × 10−16
6.7 × 10−15
2.4 × 10−14NC
8.9 × 10−16
2.2 × 10−15
5.2 × 10−15
7.1 × 10−15(65)

なお，Rc = L0 /10,  = 0.02, U∞ = 0.5, xc = yc =0, ρ = 1.0 とした．
移流した渦が一周して元の位置に戻ってくるまで計算
し，z 一定の面の速度分布を図 3 に示す．ここでは GP の
次数を N = 2 としている．図 3 の NC では，一様流を保持しない (a),(d),(g) は明らかに渦が崩れており，渦の移流でも大きな誤差が生じると分かる．しかし，(j) では一様流を保持しない（表 2 参照）にも関わらず，渦を保持している．この違いは面積保存則を満たすかどうかである(10) ．K < 2N となる (a),(d),(g) の 3 ケースは面積保存則を満足せず，空間メトリックを適切に計算出来ていない．一方，K ≥ 2N となる (j) のケースでは面積保存則を満たし，空間メトリックを適切に計算出来ている．(j)のケースでも体積保存則は満足しない為に一様流保持は成立しないが，今回の計算条件では体積保存則の誤差が解像度を大きく損なう程には生じなかったと考えられる．
実際，一様流保持の L1 誤差を示す表 2，3 でも，K ≥ 2Nの時（面積保存則の誤差がなくなり，空間メトリックを適切に計算出来る場合）に誤差が数桁下がる．時間積分を 1 次精度（Euler 陽解法）とした場合には，K = 2Nの時の一様流保持の誤差は K < 2N の場合と同程度となり，渦の解像度も下がる（図 4）．この事からも，(j) のケースは時間方向の誤差を含んでいる事が分かる．
SSC と CSC に関しては一様流を保持しない (b) のケー
ス以外は，解の次数による差異はあるものの，渦を保持出来ている．注目すべきは K = 4 での (k)SSC と (l)CSCの違いである．Sec.4.2.1 で述べたように，SSC はメト
リックの補間多項式が GP の次数 N で制限されてしま

Copyright c 2014 by JSFM

第 28 回数値流体力学シンポジウム
C06-4

い，K > N の時に解の空間精度を保てないと考えられる．これは空間メトリックに対する研究結果(10) からの予想であるが，今回移動変形格子を用いた場合にもメトリッ
クの空間次数と解の空間次数を一致させた (l)CSC の方が (k)SSC よりも渦を良く保持している．

(a) NC (K = 1)(b) SSC (K = 1)(c) CSC (K = 1)(d) NC (K = 2)(e) SSC (K = 2)(f) CSC (K = 2)(g) NC (K = 3)(h) SSC (K = 3)(i) CSC (K = 3)(j) NC (K = 4)(k) SSC (K = 4)(l) CSC (K = 4)

3. 3 次元波状変形格子を用いて一様流・移流渦の検証計算を行った．一様流保持・保存量保存性は事前に予想した通りの結果を示した．移流渦については，一様流を保持しない NC を用いても渦を保持するケースが生じた．これは今回の計算条件では，体積保存則に比べ面積保存則の誤差が格段に大きく，体積保存則の誤差が移流渦の解像度に殆ど影響を与えない為と考えられた（時間方向の離散化誤差は含まれている）．また，メトリックの空間次数が解の次数に一致しない SSC は，CSC に比べて渦の解像度で劣るケースがある事を示した．

(a) NC (K = 1)(b) SSC (K = 1)(c) CSC (K = 1)(d) NC (K = 2)(e) SSC (K = 2)(f) CSC (K = 2)(g) NC (K = 3)(h) SSC (K = 3)(i) CSC (K = 3)(j) NC (K = 4)(k) SSC (K = 4)(l) CSC (K = 4)

図 3: The contours of velocity distribution on vertical(z-constant) plane are shown. The order of the shapefunction is N = 2.

結論
3 次元移動変形格子を保存型 FR 法で用いる際に，対称保存型メトリックを用いて一様流保持・保存量保存性を両立する手法を示した．

6.

1. 再構築した流束が一様流保持（体積保存則を満足）する為には，SP 及び correction part（式 (26)）がそ
れぞれ零となれば十分である事を示した．また，保存量保存性は各セル境界でメトリックの値が一意に定まる事が十分である事を示した．
2. 従来から用いられている NC（非保存型メトリック）
は，解の補間多項式の次数を上げても一様流を保持しない事が分かった6 ．これに対し，一様流を保持す
る対称保存型メトリックの構築方法を詳細に説明し，SSC と CSC の二つを導入した．
6
セル形状の時間変化を表現する形状関数の次数と，支配方程式の時間微分を表す多項式の近似次数を独立に定めれば可能になる．

図 4: The contours of velocity distribution on vertical(z-constant) plane are shown. The order of the shapefunction is N = 2. The Euler explicit scheme is employed for the time integration (ﬁrst-order scheme intime discretization).

今後の課題として，時間方向の精度を検証する事が挙げられる．また，本稿で提案した CSC は時間方向にも精度を保証するよう構築しているが，厳密な時間メトリックは
時間方向に 3T − 1 次の多項式となる．これを表現出来る
ように時間方向にも CGP を追加する手法や，3T −1 → Tの低次近似時にエイリアシング誤差除去を行う手法7 等が考えられ，解像度が向上すると期待出来る．
7

8

3N → K の低次近似を行う空間方向にも適用可能である

Copyright c 2014 by JSFM

第 28 回数値流体力学シンポジウム
C06-4

謝辞

形格子を用いた場合の Euler 陽解法は以下のように表現出来る．

本研究は JSPS 科研費 258793 の助成を受けた．ここに記して謝意を表する．

7.

J(τn+1 )Q(τn+1 ) =J(τn )Q(τn ) − ∆τ RHS(r(τn ), r(τn+1 ), r(τn ), Q(τn ))補遺 A:剛体運動をする場合の一様流保持性

Sec.4.1 で述べたように，一般の移動変形格子では NCを用いて一様流保持出来ない．しかし，セル形状（ヤコビアンの値）が時間変化しない剛体運動格子においては，空間方向の次数が十分に高ければ（K ≥ 3N の場合には）体積保存則を満足し，一様流を保持すると予想される．表
5 には，N = 1 とした場合に 3 次元波状格子(10) を剛体回転させた場合の一様流の L∞ 誤差を示す．この時，K ≥ 3
で NC は丸め誤差で一様流を保持し，本稿の定式化（式(31),(32)）通りの結果を得ている．

(68)

本稿で用いた 2 段 2 次の Runge-Kutta 法は以下の通りである．
J(τn+1/2 )Q(τn+1/2 ) =∆τ
RHS(r(τn ), r(τn+1 ), r(τn ), Q(τn )),2
J(τn+1 )Q(τn+1 ) =J(τn )Q(τn ) −(69)

J(τn )Q(τn ) − ∆τ RHS(r(τn ), r(τn+1 ), r(τn+1/2 ), Q(τn+1/2 )),
(70)

参考文献
表 5:
L∞ norm of a freestream error(max{|v − v0 |/u0 , |w − w0 |/u0 }). The order of ashape fuction is N = 1. The wavy grid is rigidlly
rotated with ω = 0.1.
Order of SP: K
1
2
3
4

8.

NC
2.1 × 10−3
3.9 × 10−5
1.9 × 10−16
6.1 × 10−16

SSC
5.6 × 10−17
1.1 × 10−16
1.5 × 10−16
4.4 × 10−16(1) Z. J. Wang, K. Fidkowski, R. Abgrall, F. Bassi,
D. Caraeni, A. Cary, H. Deconinck, R. Hartmann,
K. Hillewaert, H. Huynh, N. Kroll, G. May, P.-O.
Persson, B. van Leer, and M. Visbal, High-order
cfd methods: Current status and perspective,
International Journal for Numerical Methods in
Fluids 00 (2012) 1–42.

CSC
6.9 × 10−17
2.9 × 10−16
2.6 × 10−15
1.6 × 10−14(2) H. T. Huynh, A ﬂux reconstruction approach tohigh-order schemes including discontinuous
galerkin methods, in AIAA 2007-4079, 2007.

補遺 B:時間メトリックに現れる時間微分の取り扱い

本稿では時間積分に 2 段 2 次 Runge-Kutta 法を用いた．ここでは，時間メトリックを計算する手順を Sjögreen
らの導入した表記(12) を用いて説明する．まず，任意の物
理量 ψτ = dψ(τ )/dτ の時間差分近似を以下のように導入
する：

ψτ '

ψ(τ2 ) − ψ(τ1 ),∆τ∆τ = τ2 − τ1 .

(3) T. Haga and S. Kawai, Toward accurate
simulation of shockwave-turbulence interaction on
unstructured meshes: A coupling of high- order fr
and lad schemes, in AIAA 2013-3065, 2013.
(4) P. D. Thomas and C. K. Lombard, Geometricconservation law and its application to ﬂow
computations on moving grids, AIAA Journal 17(1979), no. 10 1030–1037.

(66)

これに倣い，時間メトリック ξˆt の時間方向の差分近似を
以下のように導入する（簡単のため，r(τ ) = r(ξ, η, ζ, τ )と略記する．
）：(5) M. R. Visbal and D. V. Gaitonde, On the use ofhigher-order ﬁnite-diﬀerence schemes on
curvilinear and deforming meshes, Journal ofComputational Physics 181 (2002), no. 1 155–185.

ξ̂t = −[{(rη × r)ζ · r + (r × rζ )η · r}τ
+{(rτ × r)η · r + (r × rη )τ · r}ζ
+{(rζ × r)τ · r + (r × rτ )ζ · r}η ]/6
¯
−1 h˘
'
(rη (τ2 ) × r(τ2 ))ζ · r(τ2 ) + (r(τ2 ) × rζ (τ2 ))η · r(τ2 )6∆τ˘¯− (rη (τ1 ) × r(τ1 ))ζ · r(τ1 ) + (r(τ1 ) × rζ (τ1 ))η · r(τ1 )˘+ ([r(τ2 ) − r(τ1 )] × r(τ3 ))η · r(τ3 )¯+ [r(τ2 ) × rη (τ2 ) − r(τ1 ) × rη (τ1 )] · r(τ3 ) ζ˘+ [rζ (τ2 ) × r(τ2 ) − rζ (τ1 ) × r(τ1 )] · r(τ3 )¯ i+ (r(τ3 ) × [r(τ2 ) − r(τ1 )])ζ · r(τ3 ) η(6) M. Vinokur and H. Yee, Extension of eﬃcient low
dissipation high order schemes for 3-d curvilinear
moving grids, Frontiers of Computational Fluid
Dynamics (2002) 129–164.

(67)

ここで，τ3 は τ1 , τ2 とは独立に定まる時刻であり，時間積分法により異なる．この時間メトリックと保存量 Q を用いて計算した支配方程式の右辺（∂ξ Ê + ∂η F̂ + ∂ζ Ĝ）を，RHS(r(τ1 ), r(τ2 ), r(τ3 ), Q) と書く．例えば，移動変

9(7) X. Deng, M. Mao, G. Tu, H. Liu, and H. Zhang,
Geometric conservation law and applications tohigh-order ﬁnite diﬀerence schemes withstationary grids, Journal of Computational
Physics 230 (2011), no. 4 1100–1115.
(8) Y. Abe, T. Nonomura, N. Iizuka, and K. Fujii,
Geometric interpretations and spatial symmetry
Copyright c 2014 by JSFM

第 28 回数値流体力学シンポジウム
C06-4

property of metrics in the conservative form forhigh-order ﬁnite-diﬀerence schemes on moving
and deforming grids, Journal of Computational
Physics 260 (2014) 163–203.
(9) T. Nonomura, N. Iizuka, and K. Fujii, Freestream
and vortex preservation properties of high-orderweno and wcns on curvilinear grids, Computers &Fluids 39 (2010), no. 2 197–214.
(10) Y. Abe, T. Haga, T. Nonomura, and K. Fujii,
On the freestream preservation of high-orderconservative ﬂux-reconstruction schemes, Journalof Computational Physics 281 (2015) 28–54.
(11) D. Kopriva, Metric identities and thediscontinuous spectral element method on
curvilinear meshes, Journal of Scientific
Computing 26 (2006) 301–326.
(12) B. Sjogreen, H. Yee, and M. Vinokur, On highorder ﬁnite-diﬀerence metric discretizations
satisfying gcl on moving and deforming grids,
Journal of Computational Physics 265 (2014),
no. 0 211 – 220.
(13) S. Kawai and S. Lele, Localized artiﬁcial
diﬀusivity scheme for discontinuity capturing on
curvilinear meshes, Journal of Computational
Physics 227 (2008) 9498–9526.

10

Copyright c 2014 by JSFM