のようになる. A ˆ (= E ˆ / Q ˆ ) などは流束ジャコビアンである. この際, 未知量を物理量 Q ˆ ではなく, その時間変化に対応する増分 Q ˆ (= Q ˆ n+1 Q ˆ n ) ととると収束判定に都合がよい.1 つの行列要素 A ˆ や B ˆ が方程式の数のマトリックス

日本計算工学会誌「計算工学」, 第 3 巻, 第 3 号, pp. 158-166, 平成 10 年 10 月 有限体積法の最前線ー高速気流計算法の最近の動向 宇宙科学研究所　　　　藤井　孝藏 1.　はじめに 1980 年代前半の物体適合座標系の導入，後半のいわゆる TVD 風上法の導入以降圧縮性流れの計算 法に関して大きな変化は起きていない．とはいえ， 実用問題に適用しようとして出てきた様々な問題 に対していろいろなアプローチが提案されている． 圧縮性流れの数値計算法の基本原理と考え方に ついては書籍[1-3]_{などを参照していただくとし，ここ} では狭議の計算法でなく，広い意味の計算法につ いて 90 年以降に起こった変化について考えてみ たいと思う．限られたスペースなので詳細は参考文 献を読んでいただくことにしたい．また，実用問題 として高レイノルズ数流れを扱うという前提で，定常 項については対流項の取り扱い，従って非粘性オ イラー方程式に話を限る． 2.　すでに利用されている技術　? 話の前提 2.1　時間積分と空間の離散化 　非粘性圧縮性流れの基礎方程式は一般座標系で次のように書ける．

∂

_Q

ˆ

∂

t

+

∂

_E

ˆ

∂ξ

+

∂

_F

ˆ

∂η

+

∂

_G

ˆ

∂ζ

=

0

(1)

ˆ

Q

は保存量からなる未知ベクトル，

E

ˆ

などはそれぞれの方向の流束ベクトルである．1 つの式で書いて あるが，実際は質量保存，運動量保存，そしてエネルギー保存の式で，2 次元なら 4 つ，3 次元であれば， 5 つの式となる．総てが保存則を記述したものであるから，一般的な有限体積法では計算格子の周りに コントロールボリュームを考え，そこからの流束の出入りを評価してその差が物理量の時間変化であると いう定義の仕方をする．このことは式(1)にコントロールボリューム

∆

v

= ∆ξ ⋅ ∆η ⋅ ∆ζ

をかけ，時間項の みを離散化して表示してみるとよくわかる．

( ˆ

Q

n+1

−

Q

ˆ

n

)

∆

v

= −∆

t ( ˆ

∫

E d

η

d

ζ

+

F d

ˆ

ζ

d

ξ +

G d

ˆ

ξ

d

η

)

(2) 差分では式(1)を，有限体積法では式(2)を扱うことになるが，計算法の基本原理に大差はないので，式 (1)を中心に話を進める．式(2)は物理面で積分を定義することもできる（文献 1 の p. 197）が簡単のためこ こでは計算面での積分を利用した．時間方向にこれらの方程式を積分していくが，その際に式(1)の右 辺をある時刻（時間ステップ）

n

で評価して時刻

n

から

n

+

1

への物理量変化

∆

Q

を計算する方法が陽 解法，右辺の評価に

n

+

1

での物理量の影響を入れるものが陰解法である． 局所的な線形化を施すと， 陰解法の基本形は

(I

+ ∆

t

∂

∂ξ

A

ˆ

+ ∆

t

∂

∂η

B

ˆ

+ ∆

t

∂

∂ζ

C )

ˆ

∆

Q

ˆ

n

_{= −∆}

t(

∂

E

ˆ

∂ξ

+

∂

_F

ˆ

∂η

+

∂

_G

ˆ

∂ζ

)

n (3) 　

η

ξ

　　　　図 1　物体適合座標系

のようになる．

A (

ˆ

=

∂

E /

ˆ

∂

Q )

ˆ

などは流束ジャコビアンである．この際，未知量を物理量

Q

ˆ

ではなく，その 時間変化に対応する増分

∆

Q (

ˆ

=

Q

ˆ

n+1

−

Q

ˆ

n

)

ととると収束判定に都合がよい．1 つの行列要素

A

ˆ

や

B

ˆ

が方程式の数のマトリックスであり，全体マトリックスはそれらを個々の要素としたブロック行列となってい る．そのため式(3)の左辺の反転には多大な演算量とメモリー容量が必要となる．そこでこれを近似的に 因数分解する方法がすでに考えられている． 　右辺の離散化については，80 年代中盤から広く使われるようになっているのが Roe や Osher に代表さ れる近似リーマン解法である．これらは局所的な線形化に加えて方程式を固有ベクトル系に分解するこ とで，スカラー方程式における風上法の優れた特徴をシステム方程式に活かした方法である．この方法 は流体力学的な特性方程式に直接的に対応しているため物理的イメージもつかみやすく急速にその利 用がすすんだ．流束制限関数，いわゆるリミターを導入し，非線形な手法とすることで，高次精度におい ても非物理的な振動を生じない工夫（TVD 法）がなされている． 2.2　一般座標の導入 　近年の数値計算には一般座標系，特に物体適合座標系が利用される．例えば図 1 のように物体を計 算格子の境界とする物体適合座標の場合，直交デカルトの計算座標上で得られる数値解が物体まわり に形成された計算格子の頂点やセルを代表する値として得られることになる．このような物体適合座標 系の利用によって複雑形状の取り扱いは以前に比べてかなり容易になった．実際の計算では隣り合う格 子はそれぞれ j-1, j, j+1 などのように序列を持っていて，計算機上でもメモリー空間で配列として連続， または一定のストライドを持っている．例えば 2 次元の格子ならば(j,k)といった具合である．このようなデ ータ配列を維持するためには序列を持った計算格子を物体まわりに作成する必要がある． さらに計算 の精度を確保するためには，「優れた」計算格子でなければならない．実際の数値シミュレーションの現 場ではこの計算格子の作成が数値計算自身や可視化による後処理に比べてはるかに大きな負担となり， 計算格子作成に 1 ヶ月，数値計算に 1 時間などといった状況が生まれている．このような背景から有限 要素法と同じように計算格子に序列を持たせない非構造格子法が流体計算でも利用されるようになって きた． 　計算格子が序列を持たない場合でも，コントロールボリュームを考え流束の出入りから物理量を計算 することができるから，有限体積法の考え方は有効で，多くの非構造格子法では有限要素法と並んで有 限体積法が利用されている．この場合，計算格子は三角形（四面体）メッシュが利用されることが多い． 高速気流の数値計算法は，構造格子の場合と非構造格子の場合で多少異なった傾向を持っているの で以下はそれぞれの場合について別々に記していこうと思う． 3.　非構造格子の世界 3.1　時間積分 　非構造格子の時間積分に関するここ数年の進歩は著しい．陽解法ではマルチグリッド（MG:多重格子 法）を利用した収束加速に対する研究が進んでいる．多重格子法の考え方は文献 1 の 46 ページを，ま たオイラー方程式における実際の利用に関しては文献 4 の第 1 巻 488 ページ第 2 巻 98 ページあたり を参照していただきたいが，基本原理は単純である．収束過程における空間の残差は格子に依存した 波数成分に分けて考えることができるが，このうちもっとも収束が遅いのは低周波，すなわち周期の大き な波数成分である．そこで多重格子法は実際利用する格子での計算，格子を一点おきに間引いた格子 での計算，さらに間引いた格子での計算などを繰り返し行いそれを 1 つのサイクルとする．これによって 広範囲の波数成分の誤差を一様に減少させる．

　多重格子法は一般にポアソン方程式のよう な単純な楕円型の問題では収束率は 0.1 程度 が得られる．つまり，マルチグリッドサイクル 1 回で残差が一桁落ちる．圧縮性流れ方程式の 基本となる双曲型の方程式の場合では，文献 5 などにあるように収束率は 0.75 程度が得られ ている．これは構造格子，非構造格子どちらの 場合も同じである．しかしながら，粘性を考慮 したナヴィエ・ストークス方程式では境界層付 近において格子が急激に小さくなるためマル チグリッドの効率は 1 桁以上も落ちてしまう．最 近の研究は非構造格子によるナヴィエ・ストー クス方程式計算でのマルチグリッド法の効率を いかに上げるかに集中している．構造格子， 非構造格子に関わらず最近注目されているの は Semi-Coarsening と呼ばれる方法である．通 常のマルチグリッドでは総ての方向に一度に 格子の削除を行って行くが，Semi-coarsening では格子の引き延ばしの大きい方向，すなわ ち物体から離れる方向を避けてそれ以外の方 向についてのみマルチグリッドを行うものであ る．これらについては理論的， 具体的な研究がある 5など)_． 非構造格子の場合には，格 子の形状などを判定して方 向性を決めていく directional coarsening という手法で実現 される．文献 5 から粘性計算 の収束履歴の例を図 2に示 す．semi-coarsening マルチ グリッドに更に前処理などを 施すことでかなりの収束加速 が実現できる． 　一方，陰解法の適用も盛ん で あ る ． 特 に Yoon & Jameson によって構造格子 で提案された LU-SGS 法[6]

は Meshov & Nakamura によって非構造格子に拡張され[7]_{，Sharov & Nakahashi}[8]_{などによって実用化が} 図られる一方，単純な 1 往復でなく繰り返しを行えば一般の GS 法に帰着できることから，さらなる効率的 なスキームの開発も模索されている[9]_{．緩和法の演算軽減化版である LU-SGS 法は非構造格子では次} にように書くことができる． 　　

∆Q

i *

=

Vol

i

∆t

i

+

(

σ

i

)

j n

2

∆S

ij













j

∑













−1

×

RHS

+

(A

j

)

i n

+

(

σ

j

)

i n

I

2

∆S

ij

∆Q

j *













j∈Lower(i )

∑













∆Q

i

= ∆Q

i *

₊

Vol

i

∆

t

_i

+

(

σ

_i

)

n_j

2

∆S

ij













j

∑













−1

×

(A

j

)

i n

+

(

σ

_j

)

_in

I

2

∆S

ij

∆Q

j













j∈Upper (i )

∑

(4) 分母がスカラー量となっている点が効率化に効いている．後に述べる構造格子の場合と比較していただ きたい．Lohner はこれに GMRES を加え，更なる収束加速を得ている[10]_{．図 3に彼の計算結果収束履} 歴の一例を示す．LU-SGS の収束もよいが，それに加えて GMRES により著しい収束の改善が図れる． 　 　　図 2　マルチグリッド収束履歴の例 　 　図 3　非構造格子における LU-SGS と GMRES による収束加速

2.2　空間の離散化 　非構造格子における問題点は， 粘性計算のときにどれだけ正確 に物体付近の粘性層を評価でき るかにある．一般に極度につぶ れた形の三角形や四面体を利用 した計算では精度が確保できな いと考えられている．実際にこの ような格子を利用して計算すると， 速度場には大きな問題がなくても 物体面における熱流束が振動す るなどといった結果になることが ある．一方で境界層，とくに乱流 の場合に層流底層内に計算格 子を入れ，しかも形のよい格子を 作るとするとセル境界の数，従って流束計算による演算負荷がべらぼうに大きくなってしまい現実的でな い．このような場合の対処として（１）物体付近は四角形または六面体格子を利用するハイブリッド法（２） プリズムやピラミッド型の格子の利用 11)_{（３）流束評価の方法を工夫することで構造格子の精度に近づけ} る方法，などが提案されている．（１）の方法は格子形状を構造格子と同じにすることで高アスペクト比の 格子を利用して構造格子と同様の精度を確保するものである．（２）は物体付近だけは何らかの格子構 造性を持たせる（もちろんデータの構造性は不要）ものである． 図 4に構造格子，プリズム格子，四面体 格子のセル面の比較とプリズム格子のイメージを示す．例えばプリズム型では物体から離れる方向はあ る種の構造格子となるから，境界層の評価は精度が確保でき，乱流モデル評価の点でも簡便化が図れ る．しかしながら，物体から離れる方向に成長するような格子を作らなければならず，本来非構造格子の 特徴であった格子形成の自動化を多少なりとも犠牲にすることになる． 3　構造格子の世界 3.1　時間積分

　時間積分はメモリー制限などから因子分解による方法（AF: Approximate Factorization 法）が利用され てきた．その傾向は変わっていない．因子分解には ADI を基本とするものと LU を基本とするものとがあ る．ADI を基本とするものはいわゆる Beam & Warming 法から派生しており，Diagonal 法，LU-ADI 法な どが利用されている．基本となる Beam & Warming 法の ADI 因子分解は次式のようになる．

(I

+ ∆

t

∂

∂ξ

A )(I

ˆ

+ ∆

t

∂

∂η

B )(I

ˆ

+ ∆

t

∂

∂ζ

C )

ˆ

∆

Q

ˆ

n

= −∆

t(

∂

E

ˆ

∂ξ

+

∂

_F

ˆ

∂η

+

∂

_G

ˆ

∂ζ

)

n (5) Diagonal 法は 3 つの因子それぞれを対角化によって近似的に評価する．例えば第 1 因子は

(I

+ ∆

t

∂

∂ξ

A )

ˆ

=

T

ξ

(

I

+ ∆

t

δ

ξ

Λ

ξ

)

T

ξ− 1 (6) のように近似化される．ここで

T

_ξ，

T

_ξ−1はそれぞれ右，左の固有ベクトルからなるマトリックスである．これ によってブロック 3 重対角行列の反転はマトリックス演算（あらかじめ書き下せるので四則演算）とスカラ ーの反転，つまりわり算に帰着できる．但し，演算を空間微分の外に出したことによって厳密な時間精度 は失われている．LU-ADI 法はさらにこれを方向別の LDU 近似因数分解（DDADI 分解）によって行き 帰りのわり算スイープに近似化したものである．1980 年以降圧縮性流体の数値計算にはスーパーコンピ ュータが主に利用されてきた．そのため，これら因子分解法はベクトル化を意識して記述されたものが多 い．なお Diagonal 法や LU-ADI 法に関しては文献 1 や文献 3 を参照されたい． 　最近，ベクトル型のスーパーコンピュータを利用しなくともかなりの大型計算ができるようになってきた． 　 構造格子（六面体格子） プリズム格子 非構造（四面体）格子 　　図 4　プリズム格子と格子のバリエーション

ワークステーションやパソコンで行われる計 算にはベクトル化は不要で，むしろ演算が単 純なものが好まれる．すでに非構造格子の 項で触 れ た が， 最近 利用 が高 まっ てい る LU-SGS 法 6)はまさにそうした方法の典型で ある．LU-SGS 法の基礎となる近似 LU 分解 は ADI 分解式(5)に対して

(I

+ ∆

t

δ

_ξb

A

ˆ

+

+ ∆

t

δ

_ηb

B

ˆ

+

+ ∆

t

δ

_ζb

C

ˆ

+

)(I

+ ∆

t

δ

_ξf

A

ˆ

−

+ ∆

t

δ

_ηf

B

ˆ

−

+

= −∆

t(

∂

ˆ

E

∂ξ

+ ∂

ˆ

F

∂η

+ ∂

ˆ

G

∂ζ

)

n (7) と書ける．固有値の正負，すなわち情報の伝播の方向が正のものと負のものに分けた分解である．近似 LU 分解の代わりに優対角を維持するために LU-ADI のような近似 LDU 分解を施すと，

(I

− ∆

t

∆ξ

A

ˆ

j, k,l−

+ ∆

t

δ

_ξ b

_ˆ

A

+

− ∆

t

∆η

B

ˆ

j,k ,l−

+ ∆

t

δ

_η b

_ˆ

B

+

− ∆

t

∆ζ

C

ˆ

j ,k ,l−

+ ∆

t

δ

_ζ b

_ˆ

C

+

)

(I

+ ∆

t( ˆ

A

+

−

A

ˆ

−

+

B

ˆ

+

−

B

ˆ

−

+

C

ˆ

+

−

C

ˆ

−

)

j, k, l

)

− 1

(I

+ ∆

t

∆ξ

A

ˆ

j, k,l+

+ ∆

t

δ

_ξ f

_ˆ

A

−

+ ∆

t

∆η

B

ˆ

j ,k ,l+

+ ∆

t

δ

_η f

_ˆ

B

−

+ ∆

t

∆ζ

C

ˆ

j ,k ,l+

+ ∆

t

δ

_ζ f

_ˆ

C

−

)

∆

Q

ˆ

n

= −∆

t(

∂

ˆ

E

∂ξ

+

∂

_F

ˆ

∂η

+

∂

_G

ˆ

∂ζ

)

n (8) となる．ここで

_A

+

=

A

+

σ

ξ

2

, A

−

₌

A

−

σ

ξ

2

なる近似を導入すると，マトリックスの反転は総てスカラー の反転，すなわち単純なわり算に帰着する（このような乱暴な近似化が意外に功を奏するのは驚きであ る）．総てのインデックスにおいて回帰性を持っているからこのままではベクトル化はできない．これをベ クトル化するためには図 5に示すハイパー面，すなわち

j

+

k

+

l

=

const.

なる面（ここでは 2 次元なの で

j

+

k

=

const.

の線）を利用してループを組む必要がある．風上差分を利用するから行きのスイープ では A 面の各点の物理量は互いに独立であり，A-1 面や A-2 面の物理量から決まる．同様に帰りのル ープでは A 面の物理量は A+1 面や A+2 面で決まる． 　いずれにしても LU-SGS 法は現在の陰解法ではもっとも演算量が少なく，しかも簡単にプログラミング できる．演算の実体はわり算と各要素に対するマトリックス演算（あらかじめ書き下せるので四則演算）だ けで，ブロック行列などを扱う必要は全くなく，緩和法の一種とも理解できることからわかるように陽解法 の延長として考えることができる．但し，格子幅が急激に変化したり，歪んでいたりする「悪い格子」に対 して収束性が落ちるなどの欠点も指摘されている（文献 3，大林，188 ページ）． 　ベクトル化が面倒で，並列化が絶望的といえる LU-SGS 法であるが，その考え方を方向別に利用すれ ば，現在ある ADI タイプの因子分解法のコードのプログラム構造を活かしつつ簡単に演算を軽減化し安 定性を増すことができる[12,13]_{．例えば， 方向の演算子に対して 1 次元 LU-SGS の近似化を導入すると，} その方向の演算は 　　 A A +1 A +2 A -1 A -2 j+k=const. j k 　　　図 5　LU-SGS におけるベクトル化面

(I

+ ∆

t

∂

∂ξ

A )

ˆ

=

(I

− ∆

t

∆ξ

A

ˆ

j, k,l −

_{+ ∆}

_t

_δ

ξb

A

ˆ

+

)(I

+ ∆

t( ˆ

A

+

−

A

ˆ

−

))

−1

(I

+ ∆

t

∆ξ

A

ˆ

j, k,l +

_{+ ∆}

_t

_δ

ξf

A

ˆ

−

)

(9) となる．この場合でも LU-SGS の流束ジャコビアンに対する近似化を導入すると演算はすべてスカラーに なる．単純に Beam & Warming 法や LU-ADI 法の各方向の演算を式(9)で置き換えればよい．演算の効 率化ができ，プログラム構造は変わらないから，ベクトル化，さらには並列化も可能となる． いわば ADI-SGS 法とでもいうべきこの方法の収束性は図 6 のようになっている．これは鈍頭物体を過ぎる超音速流 れについての計算結果である．また，超音速においては風上側から総ての物理量が決まるので，この考 え方をうまく使って空間前進による効率解法も形成できる(Space Marching 法）[13]_． 　上記のような緩和法的な陰解法に対して，陽解法の拡張も進んでいる．すでに述べたように構造格子 においてもマルチグリッド法の利用が進んでいる．また，さまざまな前処理法が実用に供されるようになっ てきた．いわゆる陰解法と工夫された陽解法のどちらがよいのか議論の分かれるところではあるが，陰解 法の LU-ADI や LU-SGS のように簡素化すると緩和法の一種と考えることができるため，陰解法と陽解法 の区別もあいまいなものになりつつある． 3.2　対流項の離散化 　最初に述べた近似リーマン解法は，高次精度化に工夫をこらして解の単調性を維持しており，一般に TVD 法と呼ばれている．近似リーマン解法には Roe や Osher の考え方が使われるが，例えば高次精度 において Roe 法を実現する方法として MUSCL 法や Yee による non-MUSCL 法などが知られている．こ れらの方法は一般に強い膨張流れに対して安定性が悪く，淀み点付近の強い衝撃波に対してカーバ ンクル現象と呼ばれる非物理解を呈することがある．これを解決する方法として Einfeld によるオリジナル HLLE(Harten-Lax-van Leer-Einfeld)法[14]に改良を加えた和田[15]_や大林[16]_{らの HLLEM 法が，また Liou} が AUSM 法[17]_{を提案している．AUSM 法は衝撃波のところで多少のオーバーシュートが生ずる欠点が} あったが，その後，和田らの AUSM-DV 法[18]_{や Liou 自身の AUSM+法}[19]_{などによって改善されている．} 国内でも，嶋，城之内が提案する一粒子法(SHUS 法と SFS 法) [20]_{は同じ様な特徴を持っている．ここで} は基本となる AUSM 法についてその考え方を示そう． 対流項の微分を数値流束の差分で表すとしたとき， 数値流束

F

_{1/ 2}を輸送と音波による擾乱伝播の 2 つ に分けて考える．すなわち，

F

_{1/ 2}

=

u

ρ

ρ

u

ρ

v

ρ

h

























+

0

p

0

0

























=

F

_{1 / 2}c

+

0

p

0

0

























(10) 但し，

F

1/ 2 c

₌

u

1/ 2

ρ

ρ

u

ρ

v

ρ

h

























L / R ここに記号

L / R

は

(

)

L / R

=

(

)

L

if M

≥

0

(

)

R

otherwise







10 -7 10 -6 10 -5 10 -4 10 -3 0 500 1,000 1,500 2,000 ADI-SGS LU-ADI RESIDUAL ITERATION 　　　　図 6　ADI-SGS 法の収束性 　- 3

次元鈍頭物体を過ぎる超音速流れ　-(11) という意味を持つ．輸送の分は

F

_{1/ 2}c

=

1

2

u

1/ 2

[(

ρΦ

)

L

+

(

ρΦ

)

R

]

−

1

2

u

1/ 2

[(

ρΦ

)

R

−

(

ρΦ

)

L

]

(12)

と評価する但し，

Φ =

(1, u,v,h)

Tである．すなわち Flux Vector Splitting(FVS)法と同じ形である．実際 に AUSM スキームは超音速領域においては FVS 法に帰着する．AUSM 系統の方法は流速による移流 と音波による擾乱伝播を分離して考える点に特徴がある．AUSM-DV，AUSM+などはこれを改良したも のだから後は論文を見ていたいただければ理解に問題はないであろう．SHUS は質量流束を Roe の FDS 法で評価することを除いては AUSM と同じである．AUSM 系のスキームは考え方が単純で演算量 も少なく，プログラミングも容易である．詳細はオリジナルの文献を参照していただきたい． 　Roe の近似リーマン解法では流束ジャコビアンを対角化し，スカラーの風上法を適用する．一般に流 束ジャコビアンを対角化する方法は一意でない．数学的なマトリックス変換を考えればわかるように Roe などの論文にある対角化でも陰解法で利用される Pulliam & Chaussee などの対角化（文献 1 の付録参 照）でも風上差分から出る数値粘性項の形は式の上では同じになる．しかし，高次精度化する際にリミタ ー操作をするために実際のスキームではその違いが生ずる．いわゆる Roe スキームはオイラー方程式を 固有ベクトルに分解するが，得られる個々の式は特性量に対する固有方程式となっていて物理的にリー マン不変量と対応がつく．これと異なった対角化を利用すると解をおかしくする可能性がある．近似リー マン解法では Roe の固有ベクトルが出ている文献 21 や 1 の付録にある対角化を利用することが必要で あると考えていただいた方がよい． 　スカラーの方程式から想像されるように，これら高解像度風上法は中心差分にマトリックス型の数値粘 性を加えたものと理解することもできる．このように考えて，より汎用性のある定式化を施した施した CUSP などのスキーム [22]が Jameson によって提案されている． 4.　共通技術と複雑問題へのアプローチ 4.1　低速流れへの利用 　圧縮性流れの数値計算法を低いマッハ数流れに直接適用しようとすると解の進みが極端に遅くなる． 定常問題ではいくらイタレーションしても残差が減らない．これは音速と流速の比が大きくなることから擾 乱の伝播速度が特性量によって大きく違ってくることから起こる問題である．一方，推進系の流れのよう に低速から超音速まで流れの状況が変化するような応用は多数ある．そこで圧縮性流れ解法を低速流 れにも利用しようという大きな動向がある．上記の問題は変数の取り方によって起こるもので単純に変数 を変えることで問題点を除去できるとする考え方から前処理によって各固有値に対する解の成長を一様 にする考え方など多くの論文が出てきている．矢部の CIP を利用した計算は前者の典型である．後者で は Choi & Mercle や Turkel らの研究 [23-25]_{が参考になる．計算プログラムの汎用性という観点からこれら} の研究からも目が離せない．なお，矢部の CIP 法は，基礎的な計算法として，また広範囲の流れを扱う 方法としても優れた特性を持っているが，本誌には既に矢部らの解説が掲載されている [26]_{ので本稿で} は省かせていただいた． 4.2　新たなるながれ　? 直交格子の利用 　これまで述べたどのような計算手法をとるにしても，実問題ではそれ以外にいろんな障害がある．どの ように複雑物体まわりに格子を切るか，物体が移動する問題をどう扱うか，乱流モデルは？ などである． 格子形成の自動化を目指した非構造格子にしても 3 次元の複雑物体に対して完全なる自動化は不可 能に近いと考えられている．構造格子では重合格子が広く利用されるようになってきたが，本当に複雑 な形状に対しては対応が容易ではない．そこで近年，直交格子が見直されている．どんなに物体が複 雑でもとにかく格子を直交で切ろうというものである．Berger らの論文 [27]_{に端を発した研究はその簡便さ} から研究対象となり急速に進展している．国内でも越智ら [28]_{などすでにいくつかの成果がある．直交格}

子のメリットは CAD などのデータを元に物体より 格子点が外にあるかどうかを判定しさえすれば （実はここにも落とし穴があるが．．．）物体境界 が見つかるため格子形成が不要である点にある． もう 1 つの利点は物体適合座標系で必要になっ た格子点の座標，形状メトリックスなどが不要と なり，メモリーがかなり節約できることにある．通 常，3 次元では座標として格子点の 3 倍，メトリッ クスとジャコビアンで 10 倍，あわせて格子点数 の 13 倍を座標関係のデータとして持つことにな るが，直交格子の場合はインデックスと格子幅 のみが必要でこれら総てが節約できる．簡単な 2 次元の直交格子の例を図 7に示す． 　直交格子では一般に物体に関わらず格子は 正方形（正六面体）であるから計算格子ででき るセルを物体が横切ることになる．ここでの流束， 物体境界条件をどのように与えるかが大きな課 題であるが，現在のところをカットセルと呼ばれる物体が占める体積を考慮する方法と，全く無視して近 似してしまう乱暴な方法とがある．いずれにしても物体形状が大きく変化する近傍では物体適合などに よって格子を次第に細分化していく方法をとらないと精度の保証がなくなる．直交格子は，そのほとんど がデータとして非構造データを利用しているが，一部には構造データとするものもあり，これからの議論 の的となろう． 　直交格子では当然ながら境界層を含んだ粘性流れを精度よく記述することはできない． 粘性流れに 拡張するには，物体付近は物体適合格子を利用するといった既存手法とのハイブリッド法が必要になる． その一例を図 8 に示す [29]_{．ここでは直交格子} に重ねて物体まわりにだけプリズム格子を利用 することで粘性を含んだ流れの精度よい計算を 可能にしている．直交格子とプリズム格子の間 ではデータ内挿を利用しており，一種の重合格 子であるといえる．同じことは構造重合格子法 でも利用でき，バックグラウンドの格子がたまた ま直交正方格子となっていると考えることもでき る． 5.　解法の選択 　以上，さまざまな観点から最近の圧縮性流体 の数値計算法を概観してきた．それでは，どの 解法を選ぶのがよいであろうか．現状では答え は 1 つではない．実際，私の研究室でも（１）重 合格子対応の構造格子の陰（陽）解法プログラ ム，（２）解適合格子を含む非構造格子のプログ ラム，（３）直交座標のプログラム，の 3 つを持っ ている．どのように使い分けるかについて私見を 述べようと思う．（３）の直交座標のプログラムは 最近開発を進めているものであるが，当然のこ とながら流れに剥離点が存在し，しかも流れの様子によってそれが大きく変化するような問題には適用 できない．直交格子を利用してナヴィエ・ストークス方程式を解くこともできるであろうが，その信頼性は 不十分といえる．この方法の利点は複雑形状の取り込みが容易であることと物体を離れた領域の格子解 像度が高いことにもあるから，実用的には形状は複雑であるが角張った物体など剥離点が固定されてい るものや，衝撃波の位置など流れのだいたいの概要を知る場合などに対して利用する．また，すでに例 　 　　図 8　直交格子と重合プリズム格子 　 　　　　図 7　2 次元直交格子の例

として示したように物体周りに重合格子を張るなど他の手法と組み合わせて利用する方法でより実用問 題に適用していくことを考えている．（２）は形状が複雑で物体適合の構造格子を切ることが絶望的な場 合や衝撃波が複雑に干渉する場合，特にそれが時間的に変化していく問題などに向いている．2 次元 の流れしかシミュレーションの対象としない場合は形状がどうであっても（２）の手法が最適かもしれない． 格子形成も自動化が可能でプログラムを道具と割り切ったブラックボックス的な使い方ができる．衝撃波 を扱う場合でも 1 次元的な衝撃波であれば（１）の方が適当である．（１）は 1980 年代からずっと利用して きた方法である．計算格子が何とか切れる場合は計算効率はよい．メモリー上にあるデータを演算が参 照するのに連続また一定ストライドのアクセスとなるからスーパーコンピュータにしてもスカラーの計算機 にしても演算は速い．重合格子を利用すれば，形状が多少複雑でも形状を構成する要素部品ごとに計 算するので，構成要素の変更などが簡単にできる．また，移動する物体に対する計算も簡単である． （１）は結局，格子がそれなりに簡単に切れる場合の主たる選択肢と思ってもらえばいいであろう．なお， 最近，計算格子を物体や現象とともに移動することで解の精度が上がることがわかってきた [30,31]_．爆風 や渦の輸送の問題で計算格子を移動させると空間の解像度以上に優れた解を得ることができる．アプリ ケーションによってはこのような利点を活かす工夫も必要であろう． 　以上のように，各手法はそれぞれに利点，欠点がある．実際にプログラムを組むときはどの手法が自分 になじむかなども開発時間に大きな影響を与えるから大切な要素である．これらを考慮した上で問題に 応じてどのアプローチが最適かを考えてプログラム開発を進める必要があろう． 　なお，最近では圧縮性流体の分野でも市販のソフトウェアが充実してきている．例えば，数値計算法で 著名なベルギーの Hirsch のグループが開発した FINE/Turbo and FINE/Aero の最新バージョンは， Hirsch 本人によれば斬新的なインターフェイスを持ちターボ機械に限ればどのような形でも 15 分位で計 算 格 子 が 切 れ ， す ぐ に 解 析 に 入 れ る と の こ と で あ る ． ま た ， 同 様 に 数 値 計 算 法 で よ く 知 ら れ た Chakravarthy が開発した CFD++はある種のグリッドレスの機能も取り入れた（註：圧縮性流れにおけるグ リッドレスの考え方は Chakravarthy が提案した）汎用性の高いソフトウェアで遠くない将来日本でも販売 されるであろう．現実には価格の問題などあろうがその動向に注意している必要がありそうである． おわりに

　Goldberg 乱流モデルで知られる Uriel Goldberg は本年 3 月都内で開かれた SST に関するワークショ ップで「計算法で効率が 2 倍速くなっても計算結果が変わるわけではない． 中心差分＋陽的な数値粘 性の計算と TVD 法の計算で結果がどれほど違うというのか．実際のところ，流れ解析で結果を大きく左 右するのは乱流モデルである．なぜ，もっと乱流モデルに注意を払わないのか」と述べている．計算法の 動向も大切ではあるが，彼の言葉を忘れてはいけない．この解説ではスペースの都合もあり乱流モデル には触れなかった．これについては文献 3 を参照されたい． 　CFD の技術はある程度成熟期を迎えているといえる．適用範囲，信頼できるデータは何かなどを心得 て使えば有力な解析の道具となったことは事実である．航空宇宙の分野では，例えば米国を眺めた場 合，実用に供されているソフトウェアは NASA の 3 つの研究所が作ったコード，Jameson によるコード，そ れにボーイングによるコードくらいである．特に NASA のコードは国内に限れば誰でも希望者が利用で きる体制になっている．スキームを開発している一部研究者は各大学にいるものの，基本となるプログラ ムの多くは NASA などから入手したもので，それを元にそれぞれの応用問題に適した修正を施している． ヨーロッパでは CFD 研究者と数学者との連携がよく，前処理や勾配法など新しいアイデアがどんどん企 業の持つプログラムに活かされていく土壌が出来上がっている．翻って，日本国内を見てみると，多くの 大学研究者がそれぞれ独自のプログラムを開発し，それも同じ計算法，同じ乱流モデルである．これで はいつまでたっても前に進むことができない．そろそろ思い切って大きな流れとなるようなプログラムを作 る時期に来ているのではないか．ソフトウェアベンダーがお仕着せで開発するようなソフトウェアではなく， みなが協力しあってディファクトスタンダードとなるようなソフトウェアを作りあげていかなければいけない． 学会や国の研究機関はこのような流れを作っていく義務を背負っていないだろうか． 　最新の圧縮性流れの数値計算法に関しては，2 つの国際会議の講演集が参考になる．1 つは AIAA （アメリカ航空宇宙学会）の CFD 会議で 2 年に 1 回夏に開かれる．1997 年の講演集はすでに発刊され ているので冊子としてでも CD-ROM としてでも購入できる．詳細は WEB サイト(http://www.aiaa.org)を見 ていただきたい．次回は 1999 年の夏にアメリカ東海岸の Norfolk で開かれる．もう 1 つは ICNMFD (International Conference on Numerical Methods in Fluid Dynamics)ですでに 30 年という歴史がある計算

法に特化した会議である．第 16 回が 1998 年 7 月にフランスのリゾート地 Arcachon で開かれたばかりで ある．こちらの講演集は毎回 Springer 社より Lecture Notes in Physics のシリーズの中で出版される．今回 のものは来年春以降になるであろう．ちなみにこの会議はもう 1 つの国際会議 ISCFD(International Symposium on CFD)と合体して ICCFD として 2000 年に京都で開催される． また，最近ワークショップ形式の小さな集まりが頻繁に行われているのでそれらの講演集も参考になる． ほとんどが大きな出版社から本となって出される． これから出版されるであろうものを文献リストに加えて おくので参考にされたい [32,33]． 文献 [1] 藤井孝藏，流体力学の数値計算法，東京大学出版会（1994） [2] 数値流体力学編集委員会編，数値流体力学シリーズ第 2 巻，圧縮性流体解析，東京大学出版会 （1995） [3] 大宮司久明他編，乱流の数値流体力学，東京大学出版会（1998）

[4] Hirsch, C.， Numerical Computation of Internal and External Flows, John Wiley and Sons（1988） [5] Marviplis, D. J.，Multigrid Strategies for Viscous Flow Solvers on Anisotropic Unstructured Meshes, ICASE Report No. 98-6, NASA CR-1998-206910（1998）

[6] Yoon S. and Jameson A.，An LU-SSSOR Scheme for the Euler and Navier-Stokes Equations, AIAA Paper 87-0600（1987）

[7] Menshov, I. and Nakamura, Y.，An Implicit Advection Upwind Splitting Scheme for Hypersonic Air Flows in Thermomechanical Nonequilibrium, 6th Int. Symp. on CFD, No. 2, pp. 815-820（1995）

[8] Sharov, D. and Nakahashi, K.， Reordering of Three-Dimensional Hybrid Unstructured Grids for Vectorized LU-SGS Navier-Stokes Computations, AIAA Paper 97-2102 （1997）

[9] 嶋英志，構造/非構造格子 CFD のための簡単な陰解法，流体力学講演会講演集，pp. 325-329 （1997）

[10] Lohner, R.，Capabilities and Issues on Unstructured Grid-CFD for High Speed Flight-Vehicles， Proc. International CFD Workshop for Supersonic Transport Design（1998）

[11] Nakahashi, K.，Prismatic Grid Method, pp. 87-105, CFD Review 95, John Wiley and Sons（1995） [12] Fujii, K.，Simple Ideas for the Accuracy and Efficiency Improvement of Compressible Flow Simulation Methods, Proc. International CFD Workshop for Supersonic Transport Design（1998）

[13] Fujii, K.，Efficiency Improvements of Unified Implicit Relaxation/Time Integration Algorithms and Their Applications, AIAA Paper 97-2105, to appear as AIAA Journal article（1997）

[14] Einfeldt, B.，On Godunov-Type Methods for Gas Dynamics, SIAM Journal of Numerical Analysis, Vol. 25, No. 2, pp.294-318 ︵1988︶

[15] 和田安弘，近似リーマン解法について—HLLEM スキームの改良と反応流への拡張，航空宇宙技 術研究所報告 TR-1189︵1992︶

[16] Obayashi, S. and Wada, Y.，Practical Formulation of a Positively Conservative Scheme, AIAA Journal, Vol. 32, No. 5, pp. 1093-1095︵1994︶

[17] Liou, M. S. and Steffen, C. J.，A New Flux Splitting Scheme, Journal of Computational Physics, Vol. 107, No. 23︵1993︶

[18] Wada. Y. and Liou, M. S.，A Flux Splitting Scheme with High-Resolution and Robustness for Discontinuities, AIAA Paper 94-0083︵1994︶

[19] Liou, M. S.，Progress Towards an Improved CFD Method: AUSM+, AIAA Paper 95-1701︵1995︶ [20] 嶋英志，城之内忠正，一粒子法について，日本航空宇宙学会第 25 期年会講演会講演集

︵1994︶

[21] Roe, P. L. and Pike, J.，Efficient Construction and Utilisation of Approximate Riemann Solutions, Computing Methods in Applied Sciences and Engineering, VI, pp. 499-518︵1984︶

[22] Jameson, A.，Artificial Diffusion, Upwind Biasing, Limiters and their Effect on Accuracy and Multigrid Convergence in Transonic and Hypersonic Flows, AIAA Paper 93-3359(1993)

[23] Choi, Y. H. and Mercle, C. L.，The Application of Preconditioning to Viscous Flows, Journal of Computational Physics, Vol. 105, pp. 207-223︵1993︶

Equations, Journal of Computational Physics, Vol. 72, pp.277-298︵1987︶

[25] Turkel, E., Vasta, V. N. and Radespiel, R.，Preconditioning Methods for Low-Speed Flows, ICASE Report No. 96-57, NASA CR201605（1996）

[26] 矢部孝，青木尊之，固体，液体，気体の統一解法を目指す CIP 法，計算工学，Vol. 1, No. 1 （1996）

[27] Berger, M. J. and Colella, P.，Local Adaptive Mesh Rfinement for Shock Hydrodynamics, Journal of Computational Physics, Vol. 82, pp. 64-84 ︵1989︶

[28] 越智章生，中村佳朗，川添博光，SST 周りの格子の生成，第 9 回数値流体力学シンポジウム講演 論文集（1995）

[29] 寺本進，藤井孝藏，直交格子を用いた 3 次元複雑形状まわりの流れ解析，日本機械学会第 76 期全国大会，講演論文集掲載予定（1998）

[30] Fujii, K.，Lessons Learned from the Computation of Blast Wave Propagation Using Overset Moving Grids, Progess in Numerical Solutions of Partial Differential Equations, John Wiley and Sons Ltd. to be published︵1999︶

[31] 藤井孝藏，Nijs, E.，移動格子の解像度に関する一考察，第 10 回数値流体力学シンポジウム講 演論文集, pp. 232-233, （1996）

[32] Caughey, D. A. and Hafez, M. M. Ed.，Computing in Futute II, John Wiley and Sons Ltd. to be published︵1999︶

[33] Hafez, M. M.，Progess in Numerical Solutions of Partial Differential Equations, John Wiley and Sons Ltd. to be published︵1999︶