第 26 回数値流体力学シンポジウム
講演番号 E05-4

非構造格子を用いたロバストで高精度な汎用圧縮性解析ソフトウェア開発の試み
Development of a Robust and Accurate General-Purpose Compressible Flow Solverfor Unstructured Meshes○ 中島 吉隆, (株)ソフトウェアクレイドル, 大阪市北区梅田 3-4-5, E-mail: nakashima@cradle.co.jp
吉田 一, (株)ソフトウェアクレイドル, 大阪市北区梅田 3-4-5, E-mail: yoshida@cradle.co.jp
渡邊 則彦, (株)ソフトウェアクレイドル, 大阪市北区梅田 3-4-5, E-mail: watanabe@cradle.co.jp
西川 裕章, National Institute of Aerospace, Hampton, VA, 23666, USA, E-mail: hiro@nianet.org
Yoshitaka NAKASHIMA, Software Cradle Co., Ltd., 3-4-5, Umeda Kita-ku, Osaka
Hajime YOSHIDA,Software Cradle Co., Ltd., 3-4-5, Umeda Kita-ku, Osaka
Norihiko WATANABE, Software Cradle Co., Ltd., 3-4-5, Umeda Kita-ku, Osaka
Hiroaki NISHIKAWA,National Institute of Aerospace, Hampton, VA, 23666, USA
We have added a robust and accurate compressible flow solver to general-purpose unstructured mesh CFD software, byusing recently proposed numerical methods. The density-based approach is applied for that compressible flow solver.
Rotated-RHLL method is used for convective fluxes, whereas a damping term is applied for diffusion fluxes. The
damping should be increased on the highly-skewed mesh, to keep robustness with unstructured hybrid mesh system.
Although simulations for high Mach number flows, especially involving shock wave phenomenon, sometimes become
unstable and low accuracy in the previous compressible flow solver, we demonstrate the accurate computation forshock wave analysis by using the present solver. We have confirmed that numerical results of the transonic flow aroundan airfoil can be improved from the previous solver.
１．諸言”密度ベースソルバー(Density-based solver)”や“圧力ベースソルバー(Pressure-based solver)”という言葉を用いて，数値流体解析ソルバーが二種類に大別されることがある．両者は，基礎方程式を解くにあたり質量保存をどのように満たしていくかという解き方の手順に違いがあり，ソルバーとしての性質が異なる．
圧力ベースソルバーとは，SIMPLE 法に代表される圧力補正法を用いた手法で，各保存式は分離して扱われる．この手法は，市販の汎用流体解析ソフトウェアの多くで用いられている手法であり，非圧縮性流体の解析や，圧縮性流体でも比較的低マッハ数な流れ場の解析に向いた手法である．一方で，密度ベースソルバー
では，各保存式は分離せずに連立して扱われる．圧力ベースソルバーが得意とする遅い流れの解析には効率上あまり向かず，衝撃波を伴うような高マッハ数な圧縮性流れの解析に向いた手法である．
著者らが販売しているソフトウェアでは，前述に違わず圧力ベースソルバーが採用されてきたが，遷音速から超音速といった比較的高速な流れ場において，計算の不安定性や精度にしばしば問題が生じることがある．このような高速な流れ場は，航空宇宙や重機械の分野において重要であり，よりロバストで高精度な圧縮性ソルバーの開発が期待される．密度ベースソルバーを用いても高マッハ数流れにおける不安定性は難しい問題であるが，近年では高精度でかつ安定性に優れた移流スキームが提案されている(1)．
また，汎用的な流体解析ソルバーにおいては，形状品質の悪いメッシュでも安定に計算できることは運用信頼性の観点から重要である．メッシュ品質への対策として，やはり近年，拡散スキームに対して歪んだメッシュでも安定で精度の良い手法が Nishikawa
により提案されている(2-4)．
そこで，著者らは最近の高精度でかつ安定性の優れた手法を取り入れ，汎用ソルバーとして信頼性の高い密度ベースソルバーの開発を試み，市販の非構造格子系流体解析ソフトウェア SCRUY/Tetra に導入した．本報告では，従来の圧力ベースソルバーとの比較を交えながら，以下に述べる最近の手法の有効性を示す．

２．数値計算手法
SCRYU/Tetra は，四面体，五面体，及び六面体からなる非構造格子を用いた流体解析ソフトウェアであり，ソルバーのほかに自動格子生成プログラムと可視化プログラムを備える．ソルバーでは，節点中心型の有限体積法を採用している．
圧縮性流体解析では，基礎方程式に質量保存式，運動量保存式，エネルギー保存式のほか，必要に応じて乱流エネルギー式などを利用する．本報告で比較に用いる圧力ベースソルバーでは，移流項の離散化には 2 次精度の風上差分法(MUSCL 法)(5)，圧力補正には SIMPLEC 法(6)，時間進行には陰解法を用いている．密度ベースソルバーで採用した手法は，以下に示す通りである．
２．１．移流項
密度ベースソルバーでの移流項の計算には，近似 Riemann 解法として Rotated-Roe-HLL(Rotated-RHLL)法(1) を用いる．RotatedRHLL 法は，Rotated-Riemann ソルバーの一種で，Roe 法(7)と HLL法(8)を組み合わせた手法である．Roe 法は高解像度な手法であるが，その一方でカーバンクル現象と呼ばれる不安定性の問題がある．一方，HLL 法はロバストな手法であるが，粘性計算の特に境界層の計算で解が散逸する傾向がある．
Rotated-RHLL 法では，Roe法と HLL 法を組み合わせることで，カーバンクル現象を抑制し，さらに境界層での解の散逸を抑える解法を実現している．
2 次精度とするため，各節点での勾配を用いて線形解を構築する．また，不連続解の近傍などで起こり得る解の振動を抑えるた
め，勾配の評価には Venkatakrishnan の制限関数(9)を適用する．
２．２．拡散項
拡散項の計算には，Nishikawa により提案された手法(2-4)を用いる．拡散項は，移流項と同じく隣り合う離散点間の流束(拡散フラックス)として与えることで評価する．本手法では，この拡散フラックスに含まれる解の勾配の算出において，離散点で評価された勾配の算術平均に加えて，離散点間の中点における左右の解の差に比例するダンピング項を考慮する．
このダンピング項の効果は，メッシュの品質が悪くなるほど大きくなるように定義されており，
1

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E05-4
非構造格子系において安定な計算となることが期待される．
なお，ダンピング項に掛かる係数α(文献 4 の式(2.21)参照)は，本計算ではα=4/3 を採用している．

３．計算結果まず，衝撃波を伴う超音速流れの非粘性解析から，移流スキームを検証した結果を示す．次に，拡散スキームについてダンピング項の有効性を示す．続いて，遷音速流れについて乱流解析を検証した結果を示してから，最後に三次元問題での検証について述べる．

２．３．時間進行法
時間進行法には，二段二次精度 Runge-Kutta 法，または二次精度の二重時間刻み法(Dual-time stepping method)を用いる．
Runge-Kutta 法は，安定性と精度に優れた手法であるが，陽解法であり，安定性条件のために各サイクルでの時間間隔を大きく取ることができない．そのため，例えば定常場を得るための解析では，膨大なサイクル数が必要となる問題は明らかである．
二重時間刻み法では，通常の時間間隔(物理時間間隔)に加えて，疑似時間間隔を用いる．各サイクルにおいて，解が物理時間間隔進行後に二次精度後退差分を満足するよう，疑似時間を用いた反復計算による補正計算を行う(10)．この反復計算における疑似時間での時間積分には，二段二次精度 Runge-Kutta 法を用いる．このとき，疑似時間間隔には，指定クーラン数で各節点毎に規定される局所時間間隔を用いる(Local time-stepping)．二重時間刻み法では，物理時間進行が陰的に扱われるため，物理時間間隔そのものはクーラン数に制限されない．そのため，時間間隔をより大きくとることができて，計算効率の向上が期待される．

３．１．移流スキームの検証
Fig. 1 と Fig. 2 に，衝撃波を伴う超音速流れの解析例を示す．どちらの問題についても，拡散項の影響を無視したオイラー方程式による解析を行っており，壁面は全てフリースリップ壁として扱う．解析メッシュは，三角形の面素を奥行方向に 1 層掃引したプリズム要素で構成される．二次元格子となる面は，ランダムな三角形非構造格子である．また，どちらの解析例も非定常現象が対象で，時間進行には Runge-Kutta 法を用いている．
Fig. 1 は衝撃波回折問題の解析結果である．Fig. 1(a)は密度の分布を示し，Fig. 1(b)はシュリーレン写真を模擬するために密度勾配の大きさを描画している．解析領域左面の上半分から，マッハ数5.09 での衝撃波の入射を考える．右面と下面を流出境界として，他の境界面はフリースリップ壁とする．使用したメッシュは，約31 万要素で，二次元面上に約 16 万節点を持つ．Nishikawa andKitamura(1)は，セル中心の有限体積法により，同一条件での解析を行っている．彼らが使用したメッシュは，400×400 の直交格子で16 万節点の二次元格子である．Fig. 1(a)に示す結果を彼らの結果(文献 1 の Fig. 13(c))と比較すると，衝撃波や滑り線といった密度の不連続や，回折に伴う膨張波の様子は良く一致しており，非構造格子を用いた本計算でも，同程度の解像度のメッシュで同等な精度が得られていることがわかる．なお，Bazhenova ら(11)による実験のシュリーレン写真(文献 11 の Fig. 26(d))と Fig. 1(b)を比較しても，濃淡は良く一致する．
Fig. 2 は二重マッハ反射問題の解析結果であり，密度の分布を示している．解析領域左面の流入境界より，30 度斜面に対するマッハ数 10 での衝撃波の入射を考える．右面を流出境界として，他の境界面はフリースリップ壁とする．使用したメッシュは，約 42万要素で，二次元面上に約 21 万節点を持つ．この問題は，Woodward and Colella(12)をはじめ，移流スキームの検証例題としてよく用いられる問題である．Fig. 2 に示す計算結果を Woodwardand Colella の計算結果(文献 12 の Fig. 9)と比較すると，衝撃波面の(a)(b)Fig. 1 Shock wave diffraction, (a) density contours, (b) density gradientcontours.

Fig. 2 Double Mach reflection, density contours.

2Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E05-4
形状だけでなく，3 重点から左下に伸びる滑り線の様子や，マッハステム後方に見られる渦形成についても良く一致する．
以上の Fig. 1 と Fig. 2 の結果より，本論で採用した移流スキームにより超音速流れの衝撃波問題を良い精度で安定に計算できていることがわかる．

３．３．乱流計算（１）RAE2822 翼型まわりの遷音速流れについて乱流解析を行う．主流マッハ数は 0.729 とし，迎角は 2.31 度，レイノルズ数は 6.5×106である．ここでは，圧力ベースソルバーでも解析を行い，実験結果(13)との比較から密度ベースソルバーの有効性を示す．

３．２．拡散スキームの検証
2.2 節の通り，本数値計算手法では拡散スキームにてダンピング項を用いている．このダンピング項の効果は，メッシュ形状の品質が悪くなるほど大きくなる．そこで，極端に歪んだメッシュを用いた二次元角柱まわりの流れ解析を行い，ダンピング項の有効性を示す．ここで，レイノルズ数は 200 として，マッハ数は 0.3とする．時間進行法には，二段二次精度 Runge-Kutta 法を用いる．
計算で使用したメッシュを Fig. 3 に示す．図中拡大図では，壁面近傍を Y 方向に 5 倍拡大して示している．メッシュは，角柱近傍で壁面から法線方向に密となるよう歪ませており，二次元面上の三角形の各面素のアスペクト比は約 50 である．
角柱近傍のその他の方向や遠方では非常に粗いメッシュとなっている．Fig. 4 には，ダンピングの係数をα=0, 4/3 として計算した場合について，揚力の時間履歴を示している．α=0 の場合，拡散スキームにてダンピング項は考慮されず，揚力には振動が見られ，計算が不安定であることがわかる．一方で，α=4/3 としてダンピング項を考慮することで，同一のメッシュであっても計算は安定化する．
通常，Fig. 3 に示すような極端なメッシュが用いられることは稀だが，例えば複雑形状において局所的に歪んだメッシュが生成された場合に，安定な計算をするために有効な手法である．

(a)

(b)

Fig. 5 RAE2822 airfoil, pressure contours, (a) pressure-based solver, (b)density-based solver.

Fig. 3 Highly-skewed mesh. In the enlarged view, the Y-axis is stretchedby a factor of 5 to show the details of the grid.

Fig. 4 Time histories of lift force calculated by using highly-skewed
mesh.

Fig. 6 RAE2822 airfoil, pressure coefficients on the airfoil.
3Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 E05-4
本解析は定常現象が対象であり，時間進行には二重時間刻み法を用いる．圧力ベースソルバーと密度ベースソルバーの計算で用いるメッシュは同一のもので，約 4 万要素，二次元面上に約 2 万節点の比較的粗いメッシュである．乱流モデルには，SST k-ωモデル(14)を用いている．
圧力ベースソルバーと密度ベースソルバーの計算結果について，翼近傍の圧力分布を Fig. 5 に示す．両者を比較すると，大局的な分布は似ているが，翼上面の後方寄りに形成される衝撃波面の様子は異なり，密度ベースソルバーでは衝撃波面がより急峻に捉えられていることがわかる．Fig. 6 では，翼表面の圧力係数を実験結果(13)と比較して示している．密度ベースソルバーでは，翼上面での衝撃波の発生を精度良く捉えることで，衝撃波発生に伴う圧力降下を良く捉えられており，同一メッシュによる解析であっても，従来解法に比べてより実験に近い結果が得られていることがわかる．

４．結言
非構造格子系汎用流体解析ソフトウェアSCRYU/Tetraに密度ベースソルバーを導入するにあたって，近年提案された手法を用いることで，ロバストで高精度な圧縮性解析ソルバーの開発を試みた．
従来手法で安定性や精度が問題になることのあった衝撃波を伴う超音速流れについて，移流スキームの検証問題として良く扱われる問題を解析し，安定して精度の良い結果が得られることがわかった．また，本手法で採用した拡散スキームについて，歪んだメッシュにおけるダンピング項の有効性を示した．最後に，翼型まわりの遷音速乱流について従来解法との比較を行い，同一メッシュであっても衝撃波を急峻に捉えることでより精度良く流れ場を捉えられることを示し，総じて本圧縮性ソルバーで採用した手法の有効性が確認できた．
参考文献(1) Nishikawa, H. and Kitamura, K., “Very simple, carbuncle-free,boundary-layer-resolving, rotated-hybrid Riemann solvers,” J.
Comp. Phys., 227 (2008), pp. 2560-2581.
(2) Nishikawa, H., “Beyond interface gradient: A general principle forconstructing diffusion schemes,” AIAA Paper, 2010-5093 (2010).
(3) Nishikawa, H., “Robust and accurate viscous discretization via
upwind scheme-I: Basic Principle,” Comp. and Fluids, 49 (2011),
pp. 62-86.
(4) Nishikawa, H., “Two Ways to Extend Diffusion Schemes toNavier-Stokes Schemes: Gradient Formula or Upwind Flux,” AIAAPaper 2011-3044 (2011).
(5) van-Leer, B., "Towards the Ultimate Conservation Difference
Scheme. 4, A New Approach to Numerical Convection," J. Comp.
Phys., 23 (1977), pp. 276-299.
(6) Patankar, S. V., "Numerical Heat Transfer and Fluid Flow,"
Hemisphere Publishing, New York (1980).
(7) Roe, P. L., "Approximate Riemann solvers, parameter vectors anddifference schemes," J. Comp. Phys., 43 (1981), pp.357-372.
(8) Harten, A., Lax, P. D., and van-Leer, B., "On upstream differencing
and Godunov-type schemes for hyperbolic-conservation laws,"
SIAM Rev., 25 (1983), pp. 35-61.
(9) Venkatakrishnan, V., “Convergence to steady state solutions of theEuler equations on unstructured grids with limiters,” J. Comp. Phys.,
118 (1995), pp.120-130.
(10) Jameson, A., “Time dependent calculations using multigrid, withapplications to unsteady flows past airfoils and wings,” AIAA Paper,
91-1596 (1991).
(11) Bazhenova, T. V., Gvozdeva, L. G., and Nettleton, M. A., “Unsteadyinteractions of shock waves,” Prog. Aerospace Sci., 21 (1984), pp.
249-331.
(12) Woodward, P. and Colella, P., “The numerical simulation oftwo-dimensional fluid flow with strong shocks,” J. Comp. Phys., 54(1984), pp.115-173.
(13) Cook, P. H., McDonald, M. A., and Firmin, M. C. P., “Aerofoil
RAE2822 - pressure distributions, and boundary layer and wakemeasurements,” AGARD Report AR 138 (1979).
(14) Menter, F. R., “Zonal two equation k-ω turbulence models for
aerodynamic flows,” AIAA-93-2906 (1993).
(15) http://aaac.larc.nasa.gov/tsab/cfdlarc/aiaa-dpw/(16) Spalart, P. R. and Allmaras, S. R., “A one-equation turbulencemodel for aerodynamic flows,” AIAA Paper, 92-0439 (1992).

３．４．乱流計算（２）AIAA が開催する Drag Prediction Workshop (DPW)(15)は，航空機の空力特性に関する数値予測のワークショップであり，最近ではNASA Common Research Model (NASA-CRM)の航空機形状について，遷音速飛行時の空力特性に関しての調査が行われている．過去に開催された DPW では，様々なメッシュやソルバーによる数値結果が報告されており，実測との比較や他の数値結果との比較がなされている．
本数値計算手法については，3.3 節で示した通り二次元翼型まわりの遷音速流れにて，従来手法からの改善が見られている．そこで，より詳細な検証として，DPW で扱われている NASA-CRM形状についての解析を行う．Fig. 7 は，主流マッハ数 0.85，迎角 2度，レイノルズ数 5×106 として解析を行った結果であり，機体表面の圧力とマッハ数が 1 の等値面を描画している．解析には 364
万要素，321 万節点のメッシュを用い，乱流モデルにはSpalart-Allmaras モデル(16)を用いており，解析自体は安定に行えている．結果の妥当性については，過去の DPW にて報告されたデータとの比較から精査し，今後開催される DPW への参加を目指している．

Fig. 7 NASA-CRM, body surface pressure and iso-surfaces of Mach 1.

4Copyright © 2012 by JSFM