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

陰解法の時間刻みと収束率が圧縮性非定常流体解析の精度・効率に与える影響
Effects of time step and convergence rate of implicit schemeon unsteady compressible flow analyses○青野光，ISAS/JAXA, 神奈川県相模原市中央区由野台 3-1-1, aono@flab.isas.jaxa.jp
岡田浩一，（株）菱友システムズ（現：（株） IDAJ） , 愛知県名古屋市港区大江町 6-19,
okada@flab.isas.jaxa.jp
野々村拓, ISAS/JAXA, 神奈川県相模原市中央区由野台 3-1-1, nonomura@flab.isas.jaxa.jp
河合 宗司, ISAS/JAXA, 神奈川県相模原市中央区由野台 3-1-1,kawai@ flab.isas.jaxa.jp
藤井孝藏，ISAS/JAXA, 神奈川県相模原市中央区由野台 3-1-1, fujii@flab.isas.jaxa.jp
Hikaru Aono, ISAS/JAXA, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210
Koichi Okada, Ryoyu systems Co. LTD, (Currently, IDAJ Co. LTD ), 6-19 Oye-cho, Minato-Ku, Nagoya ,Aichi, 455-0024
Taku Nonomura, ISAS/JAXA, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210
Soshi Kawai, ISAS/JAXA, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210
Kozo Fujii, ISAS/JAXA, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210
The time accuracy and efficiency of the unsteady compressible flow analyses using implicit time integration are
discussed. The 2nd order implicit backward difference formulation (I2BDF) and 2nd order optimized implicit backward
difference formulation (I2optBDF) are compared. In both I2BDF and I2optBDF, 3 or 4 digit convergence of relative
residual and max CFL number from 2 to5 are best combination in viewpoint temporal order and computational times.
The difference of I2BDF and I2BDF are a little.
１．序論
近年高解像度スキームを用いた，圧縮性の音響や渦の解析が行われるようになり，その有効性が示されてきている．特に壁がある流れ場の場合は，厳しい CFL 条件により時間積分法には陽解法
よりも，陰解法が使われることが多いが，その陰的時間積分法自体の時間精度については，空間精度に比べると今まであまり検証や議論はされてこなかった．特に陰解法は巨大な行列反転を避けるために，流束ジャコビアン行列を対角化したり因子化する過程
で様々な誤差を含み実質的には時間 1 次精度未満であるので，非定常計算では内部反復法を用いるがその回数や収束性についての検証や議論は殆ど無い．
既往研究としては Choi と Moin(1)によって，完全陰解法を用いて非圧縮のチャンネル乱流を解く上で正確な解析を行うためには，
時間積分幅を dt+<0.4 (最大 CFL=3.0：非圧縮なので一様流速と格子幅で定義された CFL 数)以下にする必要があり，大きすぎる dt
を用いると乱流変動を保てずに流れ場が層流化することを示したが，内部反復の収束性については議論していない．今井ら(2)はオイ
ラー陰解法を用いて 2 次元の円柱によるカルマン渦の放出を解いた場合に，内部反復を用いなくても，境界層付近で O(1)のオーダ
ー，すなわち渦の成長する粘性の弱い領域では局所 CFL 数(一様流音速と格子幅で定義された CFL 数，以下特に断りが無ければCFL 数はこの定義による)が 1 程度であるような，十分に小さい dt
を用いれば十分に妥当な解を得られることを示したが，収束性の影響については議論していない．
寺島ら(3)は遷音速フラッターにおいては，十分に小さい時間刻み幅（振動 1 周期を 4000 分割する程
度の時間刻み幅）を用いれば内部反復法を用いなくても，内部反復法を用いた解と変わらない結果が得られること，大きい時間刻
み幅では内部反復法が有効であることから，時間精度自体を向上させるためことよりも，時間項を含めた方程式への収束性を高めることの方が重要であることを示したが，一番小さい dt の場合でも最大 CFL 数が 40 と大きく，対象も我々の興味がある音響や渦
流れではない．Bijl ら(4)は層流の 2 次元円柱流れに対して，圧縮性
ソルバーを用いて各種陰的陽解法(EDSIRK)と陰解法の比較を行
Copyright © 2014 by JSFM1

い，4 次精度の EDSIRK が精度と計算時間共に効率的であることを示したが，我々の興味がある乱流場へは適用されていない．
一方近年 Housman ら(5)よって圧縮性ソルバーを用いた音響波の伝播解析においては，内部反復法の回数による許容誤差が，解析結果に大きい影響を持つことを示し，内部反復回数との関係から，
大きい dt を用いることが必ずしも計算コスト減らすことに繋がる
わけでは無いことも示した．また野々村ら(6)は Housman らと同様な 2 次元音響波の解析をパラメトリックに行い，音響波の伝播においては計算が収束した解のうち最大CFL=2 で内部反復が2 回な
いし最大 CFL=5 で内部反復が 5 回の場合が，計算コストが小さく最適であると結論づけたが，解の時間精度については議論していない．
我々の解析(7)から経験的には最大 CFL 数が 1 から 2 程度で，非定常残差が 3，4 桁程度収束すれば良いことが分かっているが，それに対して dt や内部反復数を変えて定量的に調査したことはない．
そこで本研究では，計算コストが効率的で且つ正確な解析を行うために必要な許容誤差を満たすような，最大 CFL 数の及び内部反復数の基準を得るために，それらをパラメトリックに変えた解析を行った．計算対象としては，特に乱流を含むような剥離流れにおける時間積分法の精度を検証するために，研究の初期段階としてまず 2 次元の円柱剥離流れに対してして解析を行った．

２．解析モデル（１）支配方程式
支配方程式として 2 次元圧縮性 Navier-Stokes 方程式を用いた．
質量保存則，運動量保存則，エネルギ保存則からなり，一様流条
件及び円柱直径 d で無次元化を行った．3 つの基本的な無次元パラメータ Re， M∞， Pr はレイノルズ数，一様流マッハ数，プラ
ントル数であり，以下のように定義される．
𝑅𝑒 =𝜌∞ 𝑢∞ 𝑑𝜇∞𝑢

𝜇 𝑐

， 𝑀∞ = 𝑎∞ ， 𝑃𝑟 = 𝑘∞ 𝑝
∞

∞

(1)

ここで μ， d， a， c p そして k は，それぞれ粘性係数，円柱直径，

音速，定圧比熱であり，添字“∞”は，一様流状態を表している．
実際の計算では，これらを一般座標へと変換したものを用いた．

第 28 回数値流体力学シンポジウム
C08-2
た．u，v はそれぞれ x，y 方向の速度，ξx，ξy，ηx，ηy はそれぞれ各方向のメトリックスである．
𝐶𝐹𝐿 = max [|𝑢𝜉𝑥 + 𝑣𝜉𝑦 | + 𝑎√𝜉𝑥2 + 𝜉𝑦2 ， |𝑢𝜂𝑥 + 𝑣𝜂𝑦 | +（３）数値計算手法
ISAS/JAXAで開発された圧縮性流体解析ソルバーLANS3D(8)，(9)を用いた．
以下に今回用いた手法を示す．
2次元圧縮性Navier-Stokes
方程式を支配方程式とし，一様流の音速a∞と代表長である円柱直径dで無次元化した．渦を効率的に解くために，メトリックス，ジャコビアン，対流項および粘性項の離散化には6次精度のコンパク
ト差分法(10)を用いた．高周波の数値振動を抑えるために10次精度
3重対角フィルター(11)を使用した（フィルター係数は0.495）．フィルタリングは毎計算ステップごとに1回行っている．時間積分法には2次精度3点後退差分(以下I2BDF)または2次精度最適化4点後退
差分(12)(以下I2optBDF)に，内部反復(13)を入れたADI-SGS陰解法(14)を用いた．nステップとn+1ステップの間に内部反復(回数をmで表す)を導入して，陰解法を一般的に表すると以下のようになる．
(𝐼 +
−

𝜕
𝜕𝑦

∆𝑡
𝑎1

∆𝑡 𝜕
𝑎1 𝜕𝑥

𝜕𝐸 (𝑚−1)

(( )
𝜕𝑄

)+

∆𝑡 𝜕
𝑎1 𝜕𝑦

𝜕𝐹 (𝑚−1)(( )𝜕𝑄𝑎 ∆𝑄(𝑚−1) +𝑎0 ∆𝑄𝑛 +𝑎−1 ∆𝑄𝑛−1 +𝑎−2 ∆𝑄𝑛−2，( 1(𝐹 (𝑚−1 ) −

∆𝑡

1
𝑅𝑒

(

𝜕

𝜕𝑥

(𝐸𝑣𝑖𝑠 (𝑚−1) ) +

𝜕
𝜕𝑦

+

)) ∆𝑄 (𝑚−1) =𝜕𝜕𝑥(𝐸 (𝑚−1) ) +(𝐹𝑣𝑖𝑠 (𝑚−1 )))(2)

ここで∆𝑄 (𝑚−1) が 0 になったとすると，(2)式の右辺すなわち非定
常 Navier-Stokes 方程式が満たされることになる．内部反復が収束しさえすれば，左辺に関わらず完全な時間 2 次精度（定式化が 2
次精度の場合）が確保される．なお各陰解法の係数を以下の表に示す．
Table 1 Coefficients of Implicit Backward Differencing Formula(BDF) Schemes．
1次精度後退差分（オイラー陰解法）
2次精度後退差分(I2BDF)
2次精度最適化後退差分(I2optBDF)

a-2
0
0
-1/6

a-1
0
1/2
1

a0
-1
-2
-5/2

a1
1
3/2
5/3

また後に示す参照解には4段4次精度古典的Runge-Kutta法（以下RK4）を用いた．内部反復数と無次元時間刻み幅については，後述する．流入境界は一様流に固定，流出境界は静圧は一様流に固定し，その他の物理量は外挿とした．壁面は滑り無しの断熱条件を用いた．乱流モデルを用いない層流計算である． 計算は一様流条
件（インパルシブスタート）から開始した．陰解法の初期解には1
ステップ前と2ステップ前の解が必要であるので，それにはRK4の解を用いた．
（４）計算格子
計算格子は O 型格子を用い，外部境界は円柱を中心に円柱直径d に対して 1000d とした．計算格子は周方向 251 点，壁面に垂直
方向に 212 点あり総格子点数は約 5 万点程度，最小格子幅は境界層の厚み方向に 1.0×10-3 である．コンパクト差分法の空間解像度は高く，我々の経験では渦流れを例にとると，計算効率は 50 倍以上の格子点数を用いた汎用的な 2 次精度数値計算に匹敵する(15)，(16)．
（５）計算条件
一様流マッハ数は計算効率の点から，圧縮性の無視できる範囲
でできるだけ大きい値を取るために M∞=0.2 とし， 一様流レイノ
ルズ数は円柱直径を代表長として Red=1000 とした．これらの条件が今井らの研究(2)を参考にした．CFL 数は以下のように定義しCopyright © 2014 by JSFM2𝑎√𝜂𝑥2 + 𝜂𝑦2 ， ](3)

無次元時間刻み幅については以下のテーブルに示すように最大
CFL 数が 1，2，5，10，20，50，100 となる 7 ケースを選び，それを後に示す相対残差の 3 つの基準値として適用した．それぞれ
相対残差を満たして収束するのに要した内部反復数と，後に示す
参照残差も記す．またそれを更に I2BDF と I2optBDF に対して行った．また参照解として RK4 で CFL=1 の計算も行ったので計 43
ケースである．
Table 2 Computational Cases．
時間積分法
4段4次精度古典的RK法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法
2次精度最適化陰解法

無次元時間刻み幅 最大CFL数 桁落ち 内部反復数
0.001
1
0.001
1
3
2
0.002
2
3
2
0.005
5
3
3
0.01
10
3
3
0.02
20
3
3-4
0.05
50
3
5-6
0.1
100
3
11-21
0.001
1
4
3
0.002
2
4
3
0.005
5
4
3
0.01
10
4
4
0.02
20
4
5-6
0.05
50
4
11-17
0.1
100
4
36-55
0.001
1
5
3
0.002
2
5
4
0.005
5
5
4-6
0.01
10
5
7-12
0.02
20
5
12-24
0.05
50
5
33-52
0.1
100
5
92-126
0.001
1
3
2
0.002
2
3
2
0.005
5
3
3-4
0.01
10
3
4-5
0.02
20
3
5-6
0.05
50
3
5-7
0.1
100
3
10-20
0.001
1
4
3
0.002
2
4
3
0.005
5
4
3-4
0.01
10
4
4-5
0.02
20
4
5-6
0.05
50
4
10-15
0.1
100
4
35-50
0.001
1
5
3
0.002
2
5
3-4
0.005
5
5
4-6
0.01
10
5
6-9
0.02
20
5
11-20
0.05
50
5
31-47
0.1
100
5
84-116

参照残差
1.06E-06
2.32E-06
6.20E-06
1.23E-05
2.68E-05
7.86E-05
1.59E-04
1.06E-06
2.36E-06
6.79E-06
1.54E-05
3.75E-05
9.20E-05
1.54E-04
1.06E-06
2.09E-06
6.69E-06
1.45E-05
3.31E-05
9.88E-05
1.83E-04
9.92E-07
2.00E-06
5.22E-06
1.33E-05
2.18E-05
8.90E-05
1.76E-04
1.27E-06
2.20E-06
6.48E-06
1.47E-05
3.41E-05
8.86E-05
1.29E-04
1.05E-06
2.63E-06
7.04E-06
1.44E-05
3.24E-05
8.86E-05
1.57E-04（６）解析データ取得及び処理区間
総計算無次元時間は1000であるが，データ処理区間は準定常状態になる600から1000とした．これは揚力係数の振幅が約19周期分である．
（７）相対残差
陰解法がどれだけ収束しているかを判断するためには，一般的には残差（保存変数の変化量の L2 ノルム）の変化量で判断するが，本研究では内部反復内の残差をわかりやすくするため相対残差という指標を導入した．
j, kはそれぞれξ，η方向のインデックスを，
jmax, kmax はそれぞれξ，η方向のインデックスの最大値を示す．
これはあるタイムステップにおける各内部反復内での残差 s(j，k)
を参照残差で無次元化した物である．参照残差はあるタイムステップより 1 ステップ前の保存変数 qold1(j，k)と 2 ステップ前の保存変数 qold2(j，k)の差分を利用した．陰解法は共 10 ステップに 1回，この相対残差を用いて収束判定を行い，指定した桁だけ相対残差が落ちれば，収束したと判定して内部反復を打ち切る．

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

𝑗𝑚𝑎𝑥

∑𝑗=1
𝑗𝑚𝑎𝑥

∑𝑗=1

𝑅𝑒𝑠𝑖𝑑𝑢𝑎𝑙
𝑅𝑒𝑓e𝑟𝑒𝑛𝑐𝑒 𝑟𝑒𝑠𝑖𝑑𝑢𝑎𝑙

=

2

∑𝑘𝑚𝑎𝑥
𝑘=1 (𝑠(𝑗,𝑘))∑𝑘𝑚𝑎𝑥
𝑘=1 (𝑞𝑜𝑙𝑑1(𝑗,𝑘)−𝑞𝑜𝑙𝑑2(𝑗,𝑘))(4)2３．計算結果陰解法の時間精度を検証するために，参照解として RK4 で
CFL=1 の解を使用した．その他計算手法と計算格子は陰解法と同じである．なお本研究では陰解法以外にも，比較的大きい dt がとれる陽解法である 6 段 2 次精度最適化 Runge-Kutta 法(17)や 6 段 2
次精度非線形最適化 Runge-Kutta 法(18)も検討したが，陰解法と比較して精度や計算時間があまり変わらず，その有効性を示すことは出来なかった．精度に関しては今回用いた陰解法が内部反復が収束していれば陽解法と同等の精度が出ていたからであると．また計算時間については，陰解法が比較的少ない反復数で収束したためである．
図 1 の縦軸は I2BDF2 と I2optBDF の揚力の時間履歴から求めたストロハル数を示し，横軸に最大 CFL 数を示している．プロット
はそれぞれ相対残差が 3 桁，4 桁，5 桁落ちたケースに対応している．どのケースも CFL が小さくなるほど RK4 の CFL=1 の解，0.236 に収束していく．なお，一般的に知られている円柱のストロ
ハル数 0.2 と比べると値が大きいが，これは流れの 3 次元性があるためであり，先行研究(19)でも示されている．また，実際レイノ
ルズ数を Re=300 で計算するとほぼ 0.2 になることは確認している．
相対残差が 3 桁落ちているケースは，I2BDF2 と I2optBDF 共に
CFL 数が 5 以上で誤差が急激に増加するが，相対残差が 5 桁落ちているケースは最大 CFL 数が 100 でもストロハル数はあまり変化しない．このことから最大 CFL 数をかなり大きくとっても相対残差が 5 桁落ちていれば，または最大 CFL 数が 5 以下であれば，I2BDF2 と I2optBDF 共にストロハル数に対する影響は少ないと言える．なお，I2BDF2 の相対残差が 3 桁落ちた場合と 4 桁場合のCFL5 のストロハル数及び I2optBDF2 の相対残差が 3 桁落ちた場
合と 4 桁場合の CFL5 と 10 のストロハル数がほぼ一致しているがこれはそれぞれ Table. 2 に示すように，ほぼ同じ内部反復数で収束しているためである．他の図でも同条件化で値が一致しているのは同じ理由である．なお，I2BDF と I2optBDF を比較すると相対残差が 3 桁落ちた場合は CFL5 から 20 で，4 桁場合は CFL5 から 10で，I2optBDF の方が RK4 の解に近いが，CFL が大きくなると逆に I2BDF の方が RK4 に近くなる．

図 2 の縦軸は揚力係数の時間履歴の分散を示し，横軸に最大CFL 数を示している．どのケースもストロハル数と同様に最大
CFL 数が小さくなるほど収束していく．I2BDF2 と I2optBDF 共に
相対残差が 3 桁落ちるケースで最大 CFL 数が 10 以上で揚力係数の分散が急激に小さくなるのは，剥離せん断層の成長が遅くなる
ことで揚力の変動が小さくなるからである．一方相対残差が 5 桁落ちているケースは I2BDF2 と I2optBDF 共に CFL を大きく 100
近くに取っても，揚力の分散はあまり変化しない．このことから
ストロハル数同様に CFL 数をかなり大きくとっても相対残差が 5
桁落ちていれば，または CFL 数が 5 以下であれば，I2BDF2 とI2optBDF 共に揚力の分散の対する影響は少ないといえる．
1.200

1.100

Variance of CL

𝑅𝑒𝑙𝑎𝑡𝑖𝑣𝑒 𝑟𝑒𝑠𝑖𝑑𝑢𝑎𝑙 =1.000
I2BDF:10^-3
I2BDF:10^-4

0.900

I2BDF:10^-5
I2optBDF:10^-3

0.800

I2optBDF:10^-4
I2optBDF:10^-5

0.700

RK4
0.600
1

10

100

CFL
Fig. 2 Variance of CL against max CFL number
図 3 は I2BDF2 と I2optBDF の RK4 とのストロハル数の誤差を百分率で示し，横軸は最大 CFL 数を示している．誤差は以下のように定義する．
𝐸𝑟𝑟𝑜𝑟 =

∆𝑆𝑡
𝑆𝑡𝑅𝐾4

=

𝑆𝑡𝐼𝑚𝑝𝑖𝑐𝑖𝑡−𝑆𝑡𝑅𝐾4

(5)

𝑆𝑡𝑅𝐾4

プロットはそれぞれ相対残差が 3 桁，4 桁，5 桁落ちた物のである．I2BDF2 と I2optBDF 共に最大 CFL 数が 10 以下では相対残差
が 3 桁，4 桁落ちているケースは傾きがおおよそ 2 であるので時間精度が 2 次精度出ている．このことから今回の解析の条件下という限られた範囲ではあるが，2BDF2 と I2optBDF 共に相対残差
が 3 桁，4 桁落ちているケースでは，最大 CFL 数は 10 以下にする
ことが望ましい．また，相対残差が 3 桁，4 桁落ちているケースは最大 CFL 数が 5 程度であれば，誤差が 1％程度であることも分かった．
10.00%
1

0.260

10

100

1.00%
I2BDF:10^-3
I2BDF:10^-4

0.255

Err

I2BDF:10^-5
I2optBDF:10^-3

0.250

0.10%

I2BDF:10^-3St

I2optBDF:10^-4

I2BDF:10^-4

I2optBDF:10^-5
0.245

I2BDF:10^-5
0.01%

RK4

I2optBDF:10^-3
I2optBDF:10^-4
I2optBDF:10^-5

0.240

0.00%

CFL

0.235
1

10

100

Fig. 3 Error of Strouhal number

CFL
Fig.1 Strouhal number against max CFL number
Copyright © 2014 by JSFM3

図 4 の縦軸は I2BDF2 と I2optBDF の計算時間を，
横軸は CFL 数を示している．縦軸の単位は秒である．I2BDF2 と I2optBDF 共に

第 28 回数値流体力学シンポジウム
C08-2
相対残差が 3 桁落ちたケースは CFL 数の増加と共に計算時間は減
少して CFL=20 程度で収束するが，相対残差が 4 桁，5 桁落ちたケ
ースは CFL が大きくなるほど dt が大きくなることで総積分回数が減るので計算時間は減る方向になるが，CFL が 50 を超えると誤差を収束させるために必要な内部反復数が増えて，計算時間は増大する．また相対残差の桁落ちに比例して計算時間が増えており
特に 5 桁落そうとすると計算時間が急増することも分かる．相対残差が 4 桁，5 桁落ちたケースは最大 CFL 数が 20 から 50 付近が一番計算時間が短いといえる．また相対残差が 4 桁，5 桁落ちたケ
ースはCFL=1では，陽解法よりも計算時間がかかっている．
なお、
I2BDF2 と I2optBDF を比較した場合計算時間の差は非常に小さい．
このことから今回の解析の条件下という限られた範囲ではあるが，計算時間という点では，最大 CFL 数が 5 から 50 の間で相対残差を 3 ないし 4 桁落すのが効率的であるといえる．

I2BDF:10^-3

CPU times(sec.)

I2BDF:10^-4

I2BDF:10^-5

10000

I2BDF_CFL005

I2BDF_CFL010

I2BDF_CFL100

Fig. 5 Vorticity Contour line of 10-3
red line : Solution of RK4
blue line : Solution of I2BDF．

14000
12000

I2BDF_CFL002

I2optBDF:10^-3

I2optBDF:10^-4

8000

I2optBDF:10^-5

6000

RK4

I2BDF_CFL002

4000
2000
0
1

10

100

I2BDF_CFL005

CFL
Fig. 4 CPU times
I2BDF_CFL010

図 5，6，7 は揚力最大時の各最大 CFL 数における，瞬間流れ場の渦度を示した物である．前述したように時間精度については
I2BDF と I2optBDF を比較した場合，差が小さかったので流れ場については I2BDF のみを示す．赤い等高線が参照解である RK4 を示す．相対残差が 3 桁，4 桁，5 桁落ちたケースは青い等高線でしめす．相対残差が 3 桁，4 桁，5 桁落ちたケースはいずれも最大
CFL 数が大きくなると，参照解である RK4 の解よりも円柱の後方で渦の移流速度が遅くなる上に数値振動も増える．特に相対残差
が 3 桁の落ちているケースの最大 CFL100 は，陽解法と比べると剥離せん断層の成長が遅れと渦の移流速度の遅れている．これは陰解法の位相誤差の影響であると考えられる．一方最大 CFL 数が2 から 5 であれば 3 桁収束していれば物体近傍場はほぼ陽解法に収束している．

I2BDF_CFL100

Fig.6 Vorticity Contour line of 10-4
red line : Solution of RK4
blue line : Solution of I2BDF．

I2BDF_CFL002

I2BDF_CFL005

I2BDF_CFL010

I2BDF_CFL100

Fig.7 Vorticity Contour line of 10-5
red line : Solution of RK4
blue line : Solution of I2BDF．
Copyright © 2014 by JSFM4

第 28 回数値流体力学シンポジウム
C08-2
図 8，9，10 は図 5，6，7 の全体図である．どのケースも相対残差の値にかかわらず，最大 CFL 数 10 以上になると，円柱から離れた領域で，陽解法との渦の位置の差が目立つ．このことから相対残差が 5 桁落ちたケースであっても，今回の解析の条件下という限られた範囲ではあるが，最大 CFL 数は 10 以下にするのが望ましいと考えられる．また，最大 CFL 数が大きくなると，総計算時間に占めるファイルターの回数が少なくなるので後流で渦が減衰しにくくなることも分かった．

I2BDF_CFL002

I2BDF_CFL005

I2BDF_CFL010
I2BDF_CFL002

I2BDF_CFL100
I2BDF_CFL005

Fig. 10 Vorticity Contour line of 10-5
red line : Solution of RK4
blue line : Solution of I2BDF．
I2BDF_CFL010

I2BDF_CFL100

Fig.8 Vorticity Contour line of 10-3．
red line : Solution of RK4
blue line : Solution of I2BDF．

I2BDF_CFL002

I2BDF_CFL005

I2BDF_CFL010

I2BDF_CFL100

Fig.9 Vorticity Contour line of 10-4
red line : Solution of RK4
blue line : Solution of I2BDF．

４．まとめ
陰解法の時間刻みと収束率が圧縮性非定常流体解析の精度・効率に与える影響を調べるために，2 次元円柱周りの解析をパラメ
トリックに行い以下のことが分かった．
・今回の解析の条件下では 2 次精度陰解法と 2 次精度最適化陰解法は共に最大 CFL 数が 5 以下であれば時間精度に対する影響は少ない．
・計算時間という点では，最大 CFL 数が 5 から 50 の間で相対残差を 3 ないし 4 桁落すのが効率的である．
・2 次精度陰解法と 2 次精度最適化陰解法を比較した場合，計算精度と計算時間共に大きな差は無い．
以上のことから今回の解析の条件下では 2 次元的な渦が支配的な現象では，相対残差を十分に収束させて，最大 CFL 数を小さくする方が良いが，計算効率を高めるために最大 CFL 数が 5 程
度であれば誤差を 1%程度に抑えることが出来ると言える．
今後は計算格子を細かいものと荒いケースでも解析を行い，より一般的な知見が得られるようにする予定である．さらにより一般的な渦流れに対してへの知見を得るために，3 次元的な渦流れであるチャンネル乱流流れの解析を行う予定である．

５．謝辞
本研究には JAXA スーパーコンピュータシステム(JSS)V-systemを使用した．また時間積分法の導入には（株）菱友システムズの田中耕太氏の協力を得た．ここに記して謝意を表する．

参考文献(1) Choi, H. and Moin, P., “Effects of the Computational Time Step on
Numerical Solutions of Turbulent Flow”, Journal of Computational
Physics, (1994), 113, pp. 1-4(2) 今井, 藤井, “非定常流のシミュレーションにおける陰解法の時間刻み幅について”, 日本機械学会 1999 年度年次大会講演論文集, (1999)(3) 寺島, 藤井, ”遷音速フラッターシミュレーションにおける陰解法の時間刻み幅選定基準について”,日本機械学會論文集.
B 編, (2005), 71, pp. 2856-2863.
(4) Bijl, H.; Carpenter, M. H.; Vatsa, V. N. & Kennedy, C. A.
Copyright © 2014 by JSFM5

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

(5)

(6)

(7)

(8)

(9)
(10)

(11)

(12)

(13)
(14)

(15)

(16)

(17)

(18)

(19)

Implicit Time Integration Schemes for the Unsteady Compressible
Navier-Stokes Equations: Laminar Flow
Journal of Computational Physics, 2002, 179, 313-329
Housman, J. A , Barad, M. F. et. al., “Space-Time Accuracy
Assessment of CFD Simulations for the Launch Environment”
AIAA Paper 2011-3650, 2011
Nonomura, T., Tsutsumi, S. et. al., “Impact of Temporal and Spatial
Resolution on the Aeroacoustic Waves from a Two-dimensional
Impinging Jet”, Proceedings of the ICCFD 5 - International
Conference on Computational Fluid Dynamics, Hawaii (2012),
2012
Teramoto, S., “Large-Eddy Simulation of Transitional Boundary
Layer with Impinging Shock wave,” AIAA Journal, Vol. 38, No. 11,
pp. 2103-2112, 2000.
Fujii, K., “Progress and Future Prospects of CFD in Aerospace-Wind
Tunnel and Beyond”, Progress in Aerospace Sciences, Vol. 41, No.
6, (2005), pp. 455-470.
Fujii, K., “CFD contributions to high-speed shock-related problems”,
Shock Waves, Vol. 18, No.2, (2008), pp145-154.
Lele, S.K., “Compact Finite Difference Scheme with Spectral-LikeResolution,” Journal of Computational Physics, Vol.103, pp. 16-22,
1992.
Gaitonde, D. V. and Visbal, R. M., “Pade Type Higher-Order
Boundary Filters for the Navier-Stokes Equations,” AIAA Journal,
Vol. 38, No. 11, pp. 2103-2112, 2000.
Hixon, R. and Visbal, M. R., “Comparison of High-Order Implicit
Time Marching Schemes for Unsteady Flow Calculations”, AIAA2007-4324, (2007)
Chakravarthy, S. R., “Relaxation Methods for Unfactored Implicit
Upwind Schemes,” AIAA Paper 84 - 0165.
Nishida, H. and Nonomura, T., “ADI-SGS Scheme on Ideal
Magnetohydrodynamics”, Journal of Computational Physics Vo.
228, Issue 9, (2006), pp. 3182-3188.
Arasawa, T., Fujii, K. and Miyaji, K., “High-Order Compact
Difference Scheme Applied to Double-Delta Wing Vortical Flows”,
Journal of Aircraft, Vol.41, No.4, (2004), pp. 953-957.
Kawai, S. and Fujii, K., “Analysis and Prediction of Thin-Airfoil
Stall Phenomena Using Hybrid Turbulent Methodology”, AIAA
Journal, Vol.43, No.5, (2005), pp. 953-961.
Bogey, C., Bailly, C. “A family of low dispersive and low dissipative
explicit schemes for flow and noise computations”, Journal ofComputational Physics, (2004), 194, pp. 194-214
Berland, J., Bogey, C. and Bailly, C., “Low-dissipation and lowdispersion fourth-order Runge-Kutta algorithm”, Computers &Fluids, (2006), 35, pp. 1459-1463
Mittal, R., et.al., “A versatile sharp interface immersed boundary
method for incompressible flows with complex boundaries”,
Journal of Computational Physics, Elsevier, (2008), 227, pp. 48254852

Copyright © 2014 by JSFM6