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

ソニックブーム波形推算精度向上のための非等方細分化・疎化による解適合格子法
Anisotropic Mesh Adaptation Using Refinement/Derefinement
for Improvement of Sonic Boom Prediction○ 大木裕介, 東北大(院)，仙台市青葉区荒巻字青葉 6-6-01，E-mail:oki@ad.mech.tohoku.ac.jp
佐々木大輔，東北大(工)，仙台市青葉区荒巻字青葉 6-6-01，E-mail:sasaki@ad.mech.tohoku.ac.jp
中橋和博，東北大(工)，仙台市青葉区荒巻字青葉 6-6-01，E-mail:naka@ad.mech.tohoku.ac.jp
Yusuke OKI，Tohoku University，6-6-01 Aramaki-aza-aoba，Sendai,980-8579，Japan.
Daisuke SASAKI，Tohoku University，6-6-01 Aramaki-aza-aoba，Sendai,980-8579，Japan.
Kazuhiro NAKAHASHI，Tohoku University，6-6-01 Aramaki-aza-aoba，Sendai,980-8579，Japan.
New criteria using Shock Function for refinement and derefinement method for unstructured mesh adaptation isdeveloped to efficiently capture the shock waves in flow computations around a supersonic airplane. Since Shock
Function criteria needs filter to remove numerical errors in supersonic flow, Edge-based Error Estimation is adopted
for the filter index. Derefinement is conducted to remove unnecessary nodes created by mesh refinement. Adaptive
mesh refinement coupled with new filter index and derefinement approach is applied to an axisymmetric model insupersonic flows for validation. It is proved that this approach can improve computational cost and accuracy to predictpressure signature in near field.
1. 緒言
超音速旅客機実現への課題の一つに地上へのソニックブームによる騒音低減が挙げられる．ソニックブームは超音速機周りの衝撃波が地上へ伝播する過程で統合され，地上到達時に急激な圧力変化を伴う波形になる現象である．
ソニックブームを正確に評価するためにブーム推算精度向上が求められている．機体近傍の圧力波形の伝播解析により推算を行う際，近傍場圧力波形の精度がブーム推算精度に影響を及ぼすため近傍場圧力波形の推算精度の向上が求められる．しかし，精度向上のために解析空間内に一様に密に格子を配置すると，格子数が大幅に増大し計算コストの面で実用的ではない．
本研究では，実用的な計算コストで圧力波形精度の向上を行うため衝撃波領域における解適合細分化法[1]の開発
および検証を行う．衝撃波領域の判定には衝撃波の波面の同定に用いられる Shock Function(衝撃波関数)を指標として用いる．本指標は衝撃波が減衰する粗い格子においても衝撃波領域の判定が可能である[2]が，超音速域で用いる場
合には数値誤差の除去フィルタが必要である．このフィルタ指標を改善することによる，判定可能領域の変化及び解析全体の計算コストへの影響を検討する．また，遠方での圧力波形推算には減衰を抑えるために高い格子密度が要求され，細分化判定の誤判定による格子点数の増加量が無視できない．疎化機能を加え，過剰な細分化に対して格子点数の最適化を行うことによる格子点数への影響を検討する．

て四面体から構成される．四面体においてこのアルゴリズ
ムを適用すると，四面体計算格子の 6 辺の中で分割する辺(分割辺)を選択し，辺の中点に新たな節点を追加する．この節点と分割辺に対してねじれの位置にある辺を通る平面で，その四面体を二分割する(図 1)．二分割を行うため，細分化
領域は細分化を指示した計算格子からある程度拡散する特徴を持つ．

図 1 四面体計算格子単体の分割

2. 数値計算手法
流体解析には将来的な実機複雑形状への対応の観点から，
当研究室で開発された三次元非構造格子流体ソルバー
TAS-code[3] (Tohoku university Aerodynamic Simulation code)を用いて行う．本研究での支配方程式は三次元圧縮性 Euler
方程式で，セル節点有限体積法によって離散化される．流束評価は HLLEW 法，時間積分は LU-SGS 陰解法を用いる．
3.細分化及び疎化手法
3.1 細分化アルゴリズム
本研究で用いられた細分化アルゴリズムは，Rivara の二分
割アルゴリズム[4,5]に基づく方法である．
本研究で用いた計算格子は三次元非粘性計算のためすべ
1

3.2 疎化手法
疎化は Rivara の Derefinement[4]を取り入れ，誤判定や必要以上の細分化が行われた場合には，不要な格子点と付随する四面体を除去し分割前の状態に復元することにより格子点数の低減を行う．疎化処理に必要な情報として，各要素の細分化レベルと分割による関連情報を細分化時に付与する．
・細分化レベルの設定疎化を行う上で，各格子点が何段階細分化されたかの指標を付与する．(図 2, 図 3)
-格子点
初期格子で与えられる点の細分化レベルを 1 とする．
細分化で追加される格子点は辺の分割でのみ得られる
ため，分割する辺の細分化レベルに 1 を加えた値とする．
-格子点以外の要素
各要素が保持する格子点の細分化レベルの最大値とす
る．
・分割情報の保持
疎化は不必要と判断された格子点と関連する要素の除去を行い，除去される点が生成される前の段階に復元する．
復元時には分割元となる要素をたどる必要があるために，一つの要素は，分割元となる要素番号，分割により得られる 2 つの要素番号を保持する．三角形要素での例を図 4に示す．
Copyright © 2011 by JSFM

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


 p
Vn  V  pNode RNK = k

ShockNode RNK = k +1pEdge RNK = k +1
Edge RNK = k

図 2 格子点及び辺の細分化レベルk

Node RNK = k( p  Shock )
k

k


Vn

f shock 

Vn
 Mn
a


V


V : 速度ベクトルa : 音速p : 圧力勾配ベクトル

図 5 衝撃波関数の模式図

k+1
k

k

Face RNK = k
Face RNK = k + 1
図 3 三角形要素の細分化レベル
１１

図 6 翼面上の衝撃波(fshock=1)の可視化図
衝撃波関数は圧力勾配方向の速度成分を用いるため，圧力勾配の絶対値に比べ計算格子による数値減衰の影響が少ない．この特性から圧力や密度の勾配の絶対値を指標とするときよりも広範囲の判定が可能である．
超音速飛行時の機体近傍場圧力波形の推算精度向上が目
的 で あ る た め ， 圧 縮 領 域 (fshock>0) だ け で な く 膨 張 領 域(fshock<0)も捉える必要がある．
そこで，細分化においては fshock
に対して，

分割元
2

3

②

3

生成される要素
5

4

6

4
5 6
7
7
図 4 分割情報(三角形要素)

(3.3)

膨張領域

: f shock  -1

(3.4)

E  EthrV pVn  V  p: f shock  1, f shock  0とすることにより細分化領域の抽出を行う．細分化判定は辺ごとに行い，衝撃波の検出は fshock=1 をまたぎかつ流れ方
向に∇fshock<0 を持つ辺を細分化対象とする．
超音速域において Shock Function(衝撃波関数)を用いる際，解析空間における数値誤差による影響を無視できない[7]．
軸対称形状周りの超音速流れ(M ∞ =1.41,α=0.0)における式
3.2 の衝撃波関数 fshock を，数値誤差に対する対策を行わずに算出した際の可視化図を図 7(a)に示す．数値誤差を含むために fshock が全域に分布してしまい，衝撃波領域を同定できない．そこで数値誤差を取り除くフィルタを判定条件に追加
する．

3.3 細分化指標
本研究では衝撃波領域の細分化を行うため，格子密度による衝撃波の減衰の影響を受けづらく，広範囲の衝撃波領域の検出が可能な Shock Function(衝撃波関数)[6]により衝撃波領域の検出を行う．
この手法は図 5 に示すように衝撃波が発生する位置で，衝撃波面に対して垂直な速度成分が音速に等しくなる性質を用いたものである．圧力勾配ベクトルは図のように衝撃波に対して垂直方向を持つことから，衝撃波面に対して垂 と圧力勾配ベクトル
直な速度成分 Vn は速度ベクトルp を用いて

衝撃波

(3.5)

ここでの E はフィルタ指標，Ethr はフィルタ閾値である．フィルタにより数値誤差の影響を取り除いた可視化図を図7(b)に示す．フィルタにより衝撃波関数の判定可能領域が変化するため減衰の影響を考慮した指標の選択が必要である．

(3.1)

で与えられる．これを，音速 a で無次元化することで衝撃波関数 fshock が求められる．

f shock  M n 

Vn
a

(3.2)

衝撃波面を可視化する際には fshock=1 の等値面が用いられる(図 6)．衝撃波はこの等値面を通過時に流れ方向に減速することから，fshock=1，∇fshock<0 の条件で衝撃波を同定することができる．
2

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C06-2
４

2

2

2

1

2

５

３

６

７
2

1

4

2

3

5

6

7

は3の分割によって生成
図 8 要素の階層構造(a)フィルタリングなし(b)フィルタリング有り
図 7 衝撃波関数の可視化

以下に疎化判定の手順を示す．
1.
細分化判定の実行
2.
細分化対象を構成する格子点を抽出
3.
1., 2.に含まれない格子点を疎化対象とする

3.4 フィルタ指標
フィルタ指標には Ait-Ali-Yahia と Habashi によって提案された Edge-based Error Estimation[8](以下，Edge Error)を用い
る．１次元において格子点２点を結ぶ辺上の物理量が線形補間されるとき，正確な値 g と補間値 gh から誤差 E は

Ee  g x   gh x 

(3.6)
2

E

2
1
2 d ghhedx 2 e
120

(3.7)

3

3

2

n: 細分化レベル
細分化対象
保護対象
2
疎化可能点

2
2

(3.8)

3

2

3

3
3

2

(b)保護対象の抽出

3

2

3

2

(c)疎化対象の決定(d)除去後
図 9 疎化判定手順(2 次元)

(3.9)

1

Ee   ( x2  x1 )T M (t )( x2  x1 )dt

33(a)細分化判定の実行(3.8)式を 3 次元に拡張すると

 2 qˆ
hi , j xi x j323と表すことができる．ここでの he は辺長さ，x は 2 点間の局所座標系である．指標は辺を構成する 2 格子点上の物理量の 2 次勾配と辺長さにより算出される．(3.7)式より，辺上における誤差の RMS は(3.8)式に示すように求められ，辺における誤差は辺長さの二乗と変数の二次勾配に比例する．
RMS
e

2

2

として表され，2 次と 1 次の補間の差は

1
d 2 ghEe  xhe  x  2
2
dx

3

3

と表せる．指標の特性から，任意の閾値に対して格子が粗い物体遠方では格子が細かい近傍よりも要求される勾配が小さく広範囲のフィルタリングが可能である．細分化が必
要な領域は周りよりも Edge Error の評価値が大きい領域であるため，閾値には解析領域全体の平均値を用いる． Edge
Error の算出に用いる物理量は圧力を用いることで，圧力の誤差を含む辺を重点的に細分化する．

2 次元格子に対する疎化判定の例を図 9 に示す．図示する格子は三角形に対して複数回細分化が行われたもので，格子点に示す番号は細分化レベルである．まず細分化判定を行い細分化領域の格子点を特定する(a)．次に細分化対象の格子点が存在するのに必要な格子点の保護を行う(b)．細分化
レベルが n の格子点は n-1 の格子点を持つ要素を分割して得られるため，辺でつながる n-1 以下の格子点はすべて保護す
る必要がある．これらの処理後に残った格子点が疎化可能点であり，疎化対象として扱う．

3.5 疎化判定
疎化は細分化された格子に対して，不要な格子点と関連する辺や四面体，三角形要素を取り除く．精度に対して重要な要素が除去されることを防ぐために，予め細分化判定を行い細分化領域を構成する格子点の保護を行う．
細分化により得られる格子は階層構造を持つため，細分化対象に加えて，対象が存在するために必要な格子点を保護する．図 8 のように左から細分化が進み右の格子が得られるとき，青丸点を保護する場合を考える．青丸は 3 の分割によって得られるため 6,7，及び 3 の分割元である 1 を構成する格子点を保護する必要がある．

3.6 細分化手順
本解適合細分化手法による近傍場圧力波形精度向上は以下の手順により行う(図 10)．
1. 初期格子の定常解の解析
2. CFD結果をもとに疎化・細分化の実行
3. 再度定常解の解析
以後，格子点数の増加が認められなくなるまで 2 と 3 を 1セットとして繰り返す．細分化の収束の判定は，細分化による格子点増加数が全格子点数の 1%未満となる時とする．
細分化により追加される格子点の物理量は，分割された

0(3.10)

3

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C06-2
辺の両端の接点における物理量の平均値を与える．細分化
後の格子に細分化前の計算結果を内挿することにより収束を速めている．また，全体の細分化-流体解析の反復回数の低減のために細分化後の格子に対して，流体解析を行わずに再び細分化する処理を取り入れる．細分化処理を内部反復することにより細分化を促進する．内部反復を行う場
合，追加点周りの判定精度が低下し必要以上の細分化が行われる傾向にあるが，格子点数の増加よりも収束性の向上
により全体の解析時間が短縮できると判断した．また，疎化処理は細分化判定精度が低下する細分化処理前に行う．
初期格子にて定常解の解析

(a)細分化領域全域

疎化処理
・細分化判定・疎化判定
疎化対象の除去(b)モデル部拡大図
図 11 初期格子(対称面)
細分化処理
細分化判定

細分化

内部反復

格子点増加数が総格子点数の1%以上定常解の解析

終了

図 10 解適合細分化手順

4.

図 12 初期格子における⊿p 分布(対称面)

検証計算

本検証では，解適合格子細分化による圧力波形精度向上の検証のために，風洞実験による実験値の存在する軸対称形状の解析を行う．物体直下 h/L=5.0 における圧力波形を細分化前後で実験値と比較することにより本細分化の有効性の検証を行う．

解析に用いた初期格子の対称面を図 11 に示す．モデル
表面格子の粗さによる衝撃波の発生を防ぐために，先端部及びスティング部(0≦x≦2L)は密な表面格子を用いている．
スティング部後部(x>2L)は格子点数低減のため，粗い格子としこの領域付近を細分化領域から除外することにより不要な細分化を防ぐ．空間格子は細分化を前提として比較的粗くしている．非粘性計算のため空間格子は四面体格子のみで構成し，総格子点数は約 35 万点である．
初期格子における主流マッハ数 1.41 での解析により得られた対称面の⊿p 分布を図 12 に示す．

4.1 解析対象
本解析では文献[9]にある断面積が線形変化する軸対称形状を用いた．先鋭物体の長さを L=1 とし，実験を模擬するために後部に同半径のスティング部(長さ 3L)を設置した．
実験で用いられた解析対象の半径 r 及び断面積 A は次式で示される(l=2inch)．

r

0.02



A  0.02 xx(0  x  l )(4.1)(0  x  l )

(4.2)

4

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C06-2(b)Edge-based Error Estimation
図 13 衝撃波関数の可視化(対称面)0.6Fshock[-]

1

-3

0.4

0.5
0
-2

-1

0.3
0

1

2

3

4

5

-0.5

Fshock(EE_Fil)
Fshock(PG_Fil)
Pressure gradient

-1
-1.5
-2

x[-]

0.2
0.1

0.1
0.09

1.5

0.08
1
0.07

Fshock[-]0.5

2

Pressure Gradient Magnitude[-]

2
1.5

0.5

0.06

0

-1

0.05
0

1

2

3

-0.5

-1.5
-2

(a)対称面上(h/L=0.5)

5

Fshock(EE_Fil)
Fshock(PG_Fil)
Pressure gradient

-1

0

4

6

0.04
0.03
0.02
0.01

Pressure Gradient Magnitude[-](a)圧力勾配値0x[-](b)対称面上(h/L=2.0)
図 14 衝撃波関数と圧力勾配の比較(a)圧力勾配(b)Edge-based Error Estimation
図 15 モデル近傍の衝撃波関数の可視化(対称面)4.2 数値誤差フィルタ比較検証
衝撃波関数の数値誤差除去フィルタの指標として
Edge-based Error Estimation(Edge Error)の有効性を確認するために，衝撃波関数を用いた衝撃波の可視化においてフィルタ指標として用いられる圧力勾配の絶対値との比較を行う．初期格子の定常解に対して，それぞれの指標によりフィルタ処理を行った衝撃波関数を可視化することにより，判定可能領域の比較を行う．各フィルタ指標の閾値はEdge Error に対して解析領域全体の EdgeError の平均を，圧力勾配に対しては 0.01 を設定した．
対称面上におけるフィルタ処理後の衝撃波関数の可視化図を図 13 に示し，h/L=0.5, 2.0 における衝撃波関数値と圧力勾配を図 14 に示す．モデル近傍(h/L=0.5)においては減衰の影響が小さく圧力勾配も閾値に対して十分得られており，EdgeError(EE_Fil)と圧力勾配(PG_Fil)のどちらも衝撃波領域を鮮明に捉えることができている． h/L=2.0 では減衰の影響が現れており，圧力勾配を指標とする場合において減衰の影響で衝撃波領域が大きく制限されている．一方でEdgeError を用いた場合には，減衰の影響は小さくより遠方も含め十分に衝撃波領域を捉えることができている．この結果より，Edge Error をフィルタ指標とすることで圧力勾配に比べより遠方まで判定可能である．
物体先端近傍での衝撃波領域の可視化図を比較すると，圧力勾配では衝撃波領域として可視化されている部分がEdge Error では誤差として除去されている(図 15)．物体先端付近は特に格子が細かく，周りに比べて大きな圧力勾配が生じていたとしても，Edge Error は格子全体の平均以下となり細分化対象から除去される．この特性により圧力の 2 次勾配に対して、必要以上の細分化が抑制されると期待できる．

5

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C06-2
4.3 解析条件
細分化に用いた設定値を表 1 に示す．本解適合格子細分
化では，フィルタリングを施した衝撃波関数を指標として細分化の判定を行う．細分化領域内で指標に対する条件を満たす範囲とし，ニ分割アルゴリズムにおける分割辺の決定には細分化に伴う格子の歪みを抑制するために四面体を構成する辺において最長のものを選択することとした．
最長辺を分割辺に選択することにより，四面体表面の最大内角を分割し，細分化による格子の歪みを抑制する．細分化後の格子密度を制限する閾値として，分割限度を定める．
分割辺の長さが分割限度以上の時に分割を行うものとす
る．
圧力波形取得は実験値の存在するモデル直下 h/L=5.0 の他に，h/L=2.0 の位置でも行いモデルからの距離に応じた圧力波形の変化を確認する．細分化領域は対称面付近のモデル下方に直方体の領域を設定する(図 16)．
衝 撃 波 関 数 の フ ィ ル タ 指 標 を 圧 力 勾 配 と Edge-based
Error Estimation とで比較し，初期格子におけるフィルタ性能の確認と共に細分化に対する影響を調べる．また，疎化処理の有無により格子点数の変化及び圧力波形への影響の確認を併せて行う

除外領域

h/L=5.0

細分化領域

図 16 細分化領域の設定

表 1 細分化設定値

Case1
R3C0_PG
フィルタ指標フィルタ閾値E thr
細分化内部反復数
細分化判定
衝撃波
膨張領域
分割限度
疎化

Case2
Case3
R3C0_EE
R3C1_EE
Edge Error

圧力勾配

解析領域の平均値
30.01

f shock= 1，∇f shock <0
f shock ≦-1
0.05

無

有

0≦x≦10-6≦y≦00≦z≦1.02≦x≦10-0.5≦y≦0.50≦z≦0.5

細分化・疎化領域

判定除外領域

表 2 細分化結果

細分化・再計算回数
格子点数

初期格子
Init
約35万点

Case1
Div5e-2R3C0PG
14
約65万点

6

Case2
Div5e-2R3C0EE
10
約55万点

Case3
Div5e-2R3C1EE
8
約51万点

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C06-2(a)Case1 圧力勾配フィルタ(疎化無し) (b)Case2 Edge Error フィルタ(疎化無し) (c)Case3 Edge Error フィルタ(疎化有り)
図 17 細分化終了時の格子及び Δp 分布(対称面)3p  h  4
 
p L

3

0.03

p  h  4 0.03
 
p L

Init

Case1:Div5e-2R3C0PG

0.02

Case2:Div5e-2R3C0EE

Case1:Div5e-2R3C0PG
Case2:Div5e-2R3C0EE

Case3:5e-2R3C1EE

0.01

Init

0.02

Case3:5e-2R3C1EE

0.01

Exp
0
-1.5

-1

-0.5

0
0

-0.01

0.5

1
x  h 
 
L L



-1

-0.5

0

1
4

0.5

1


-0.01

1

x  h  4
 
L L

-0.02

-0.02

(a)対称面 h/L=2.0(b)対称面 h/L=5.0
図 18 圧力波形の比較

4.4 解析結果及び考察
初期格子に細分化を行うことで，h/L=5.0 における圧力波の伝播が改善していることが確認できる(図 17，図 18)．
しかし，実験値と比較すると波形のピーク値が減衰していることに加え，ピーク位置も後方にずれている．細分化を実行することで，衝撃波による初期圧力上昇が改善されると共にピーク位置が前方に遷移していることから，より分割限度を小さくとり格子密度を高めることで波形の推算値が改善すると考えられる．

細分化過程の違いを確認するため反復数 5 回の時の格子及びΔp 分布を図 19 に示すと共に，各反復回数での総格子点数と増加数のグラフを図 20 に示す．
圧力勾配を指標とするとき，細分化による圧力伝播の改善により遠方へ細分化領域が広がるために，反復数 5 回の段階では h/L=4.0 付近で細分化が止まっている(図 19)．一方，Edge Error を指標として用いた場合，粗い格子での減衰の影響を受けず，この段階で既に h/L=6.0(設定領域の最遠方)まで細分化が進んでいる．フィルタ指標と衝撃波関数の双方で減衰の影響を受けにくい特性を持つことから，細分化初期の段階から近傍から遠方まで均一に細分化が行われたものと考えられる．格子点数の変化履歴(図 20)から
Edge Error を用いた方が圧力勾配を用いたときよりも，格子点の増加量が多く早い段階で収束している．また，細分化終了時の格子点数が EdgeError を指標としたときに約 10
万点減少している．EdgeError の圧力変動に応じた細分化
により，圧力変動が小さい領域で強い細分化が行われない特性によるものと考えられる．
精度に対して細分化に要する反復回数と格子点数を減らせることからフィルタ指標として Edge Error が適していると考えられる．

4.4.1 数値誤差フィルタの影響
表 2 より細分化完了までに要した細分化-流体解析の反復回数は圧力勾配を用いた際には 14 回であったのに対して，Edge Error を用いた際には 10 回と減少していることが確認できる．図 17 の格子の可視化図を比較すると，フィルタ指標に圧力勾配を用いたものは近傍から遠方まで同程度の細分化が行われているのに対して，EdgeError を用いたものは遠方で近傍に比べて細分化が弱く粗い格子となっている．これは，フィルタ判定に圧力変動と格子の大きさを用いるため，分割限度の不足から遠方での圧力変動が減衰し圧力変動に応じた細分化が行われたためと考えられる．また，圧力波形の変化は h/L=2.0 においては違いがほとんど見られないが，h/L=5.0 では，フィルタ指標により若干の違いが見られる．

7

Copyright © 2011 by JSFM

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

140,000

700,000

140,000

600,000

120,000

600,000

120,000

500,000

100,000

500,000

80,000

N_P_F(Case2:R3C0EE)
N_ADD(Case1:R3C0PG)

300,000

60,000

N_ADD(Case2:R3C0EE)
200,000

40,000100,000

20,0000

総格子点数[-]

N_P_F(Case1:R3C0PG)
400,000

510
細分化回数[-]

15

80,000

N_P_F(Case3:R3C1EE)
300,000

60,000

N_ADD(Case2:R3C0EE)
N_ADD(Case3:R3C1EE)

200,000

40,000100,00000

100,000
N_P_F(Case2:R3C0EE)

400,000

20,000

0

20

0

0

図 20 格子点数履歴(フィルタ指標)

格子点増加数[-]

700,000

格子点増加数[-]

総格子点数[-](a)Case1 圧力勾配フィルタ(疎化無し)(b)Case2 Edge Error フィルタ(疎化無し)
図 19 細分化 5 回終了時の格子及び Δp 分布(対称面)510
細分化回数[-]

1520

図 21 格子点数履歴(疎化の有無)

4.4.2 疎化による格子点数削減効果
細分化終了時の格子を比較すると疎化の有無で大きな変化は確認できないものの，h/L=3.0 や細分化領域末端部分において格子密度が低くなっていることが確認できる．格子点数の変化履歴を図 21 に示す．Case3 は疎化を含むため格子
点増加数を
格子点増加数=格子追加数(細分化)-格子除去数(疎化)として算出している．反復回数に対する格子点数の変化は，疎化処理を加えることで各細分化による格子点の増加数が減少しており，収束性への影響は見られない．
格子点数は疎化により 4 万点減少しており，初期格子からの格子点増加数の約 20%に当たる．疎化処理に要する解析
コストは流体解析に比べ十分に小さいため，流体解析を含めた細分化処理全体の解析コストの削減が期待できる．

0.03

3

p  h  4
 
p L

Case3:5e-2R3C1EE

Case4:Div3e-2R3C1EE

0.02

Case5:Div1e-2R3C1EE
0.01

0
-1.5

-1

-0.5

0

0.5

1


1

x  h  4
 
L L

-0.01

-0.02

(a)対称面 h/L=2.0

4.5 分割限度の検討
分割限度を 0.05 としたときの，h/L=5.0 での圧力波形は実験値に比べ減衰していた．より分割限度を小さくした場合に圧力波形が改善されるかの検討を行う．
前項までに，Edge-based Error Estimation(Edge Error)を指標とする数値誤差フィルタと疎化機能が有効であると確認した．解析条件は，4.3 節の Case3(R3C1_EE)の設定値を基準として，分割限度を 0.03 及び 0.01 と小さく与え効果の検討を行う．細分化設定及び結果を表 3 に示す．

3

p  h  4
 
p L

0.03
Case3:5e-2R3C1EE
Case4:Div3e-2R3C1EE

0.02

Case5:Div1e-2R3C1EE

Exp

0.01

0
-1

-0.5

0

0.5

1


分割限度 0.05 は 4.3 節の Case3 の結果を用いた．Case4，Case5
では細分化が多く行われるため，細分化完了までの反復回数
がそれぞれ 13，11 回と増加した．格子点数は Case3 と 4 では大きな増加は見られなかったが，Case5 において格子点数
が 10 倍以上に急増している．

1

x  h  4
 
L L

-0.01

-0.02

(b)対称面 h/L=5.0
図 22 圧力波形の比較

表 3 細分化設定・結果

分割限度
細分化・再計算回数
格子点数

Case3
Div5e-2R3C1EE
0.05
8
約88万点

Case4
Div3e-2R3C1EE
0.03
13
約99万点

Case5
Div1e-2R3C1EE
0.01
11
約1311万点

8

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C06-2(a) Case3(Div5e-2)(b) Case4(Div3e-2)
図 23 格子可視化図(対称面上)

圧力波形と対称面上における格子の可視化図をそれぞれ
図 22，図 23 に示す．分割限度を小さく取ることで，衝撃波と膨張波の発生位置の格子密度が高くなっていることが確認できる．特に Case5 では遠方での圧力変動が Case3,Case4 に比べはっきりと伝わっており，細分化も遠方まで細かく行われている．急激な格子点数の増加は，分割限度の縮小により圧力変動の伝播が改善されることにより遠方においても分割限度まで細分化が行われたためと考えられる．また，細分化の
内部反復の細分化判定精度の低下の影響により，内部反復を行わない条件に比べ格子点数が大きく出ていると考えられる．
圧力波形を比較すると，h/L=2.0 の位置では波形のピークが高まるともに，ピーク位置も前方に移動しており改善傾向が見られる．h/L=5.0 では Case4 の条件では h/L=2.0 に比べ大きな改善は見られないが，Case5 の条件に於いては波形形状に明らかな変化が生じより実験結果に近いものとなった．この結果から，分割限度をより小さくすることで，圧力波形の推算精度の改善が期待できる一方で，より格子点数が増加すると考えられる．
5．結言
近傍場圧力波形の推算精度向上のために，Shock Function
を 指 標 と す る 解 適 合 細 分 化 法 に 対 し て Edge-based Error
Estimation(Edge Error)を指標とする誤差フィルタおよび疎化機能を取り入れた．Edge Error 指標は広範囲のフィルタに適しており，衝撃波関数の一般的なフィルタ指標である圧力勾配の絶対値に比べて精度を損なわずに細分化-流体解析の反復回数と格子点数を減らすことを確認した．疎化機能の追加
では流体解析後不要な格子点を除去でき、疎化機能を行わない条件に比べ格子点増加数を約 20%低減できることを確認した．これらの結果より，誤差フィルタ指標と疎化処理を組み合わせることにより，軸対称形状の解析で最大約 53％の格子点数の低減を確認した．
また，分割限度をより小さく設定し格子密度を高めること
により波形推算精度の改善の傾向が見られたが，十分な分割限度での解析には至らず，細分化の内部反復に伴う細分化判
定精度の低下が無視できないことを確認した．
今後は，より高い格子密度までの細分化に対応するため，少ない細分化で精度の改善が期待できる異方性の強い分割法の検討を行う．

9(c) Case5(Div1e-2)

謝辞
本研究は，JAXA からの委託研究「ソニックブーム推算のための高精度近傍場圧力波形推算技術に関する研究」により行われた．
本研究の計算結果の一部は，東北大学サイバーサイエンスセンターおよび JSS(JAXA Supercomputer System)のスーパーコンピュータ SX-9 を用いて得られた．
参考文献[1] Murayama, M., Nakahashi, K. and Sawada, K., “Simulation of
Vortex Breakdown Using Adaptive Grid Refinement with
Vortex-Center Identification,” AIAA Journal, Vol. 39, No. 7, pp.
1305-1312, July 2001.
[2] 大木裕介，ソニックブーム波形推算精度向上のためのShock Functionを用いた解適合細分化法，第42回流体力学
講演会/航空宇宙数値シミュレーシン技術シンポジウム
2010 講演集，2010[3] Nakahashi, K., Ito, Y. and Togashi, F., “Some Challenges ofRealistic Flow Simulations by Unstructured Grid CFD,”
International Journal for Numerical Methods in Fluids, Vol.43,
Issue 6-7, pp.769-783, 2003.
[4] Rivara,
M. C., “Selective Refinement/Derefinement
Algorithms for Sequences of Nested Triangulations,”
International Journal of Numerical Methods in Engineering,Vol. 28, pp. 2889-2906, 1989.
[5] Sharov, D. and Fujii, K., “Three-Dimensional Adaptive
Bisection of Unstructured Grids for Transient Compressible
Flow Computations,” AIAA Paper 95-1708, 1995.
[6] Lovely, D. and Haimes, R., “Shock Detection from
Computational Fluid Dynamics Results,” AIAA Paper 99-3285,
1999.
[7] Paparone, L. and Tognaccini, R., “Computational Fluid
Dynamics-Based Drag Prediction and Decomposition,” AIAA
Journal, Vol. 41, No. 9, pp. 1647-1657, 2003.
[8] Ait-Ali-Yahia, Habashi WG. “A directionally adaptive
methodology using an edged-based error estimate on
quadrilateral grids,” International Journal for Numerical
Methods in Fluids, 1996; 23:673–690.
[9] Carlson, H., et al. “A Wind-Tunnel Investigation Of The Effect
Of Body Shape On Sonic-Boom Pressure Distributions,”
NASA TN D-3106, 1965.

Copyright © 2011 by JSFM