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

非圧縮流れによる残差切除法を用いた高精度流れ解析法の検証
Verification by Incompressible Flow of a High Accuracy Flow Analysis Method Using
the Residual Cutting Method○ 菊地一雄, JAXA, 東京都調布市深大寺東町 7-44-1, E-mail: kikuchi@chofu.jaxa.jp
西澤敏雄, JAXA, 東京都調布市深大寺東町 7-44-1, E-mail: nishizawa.toshio@jaxa.jp
関根義人, ASIRI, 東京都千代田区内神田 1-18-14, E-mail: sekine@asiri.co.jp
Kazuo KIKUCHI, JAXA, 7-44-1 Jindaiji-higashi, Chofu, Tokyo.
Toshio NISHIZAWA, JAXA, 7-44-1 Jindaiji-higashi, Chofu, Tokyo.
Yoshihito SEKINE , ASIRI, 1-18-14, Uchi-Kanda, Chiyoda, Tokyo
The goals of the present study are to develop the numerical flow analysis method based on a high-accuracy algorithm and to simulate flow
field around realistic shape like cascades of gas turbines. In this paper, the algorithm which forms the basis of the high-accuracy flow
simulation is described. In this algorithm, firstly the pressure gradient of Poisson equation is corrected in order to make the residual of theequation of continuity in the next iteration into zero. Then, new pressure fields are calculated with the corrected pressure gradient. Finally,momentum equations are solved using new pressure field. The target accuracy is realized by repeating the above procedure. The ResidualCutting Algorithm, which JAXA has developed, is adopted as matrix solver. Three cases of incompressible flow are analyzed to verify ournumerical code; Poiseuille flow , flat plate blades, and square cavity. The results are good and the applicability of our high-accuracy method
is demonstrated.

１. はじめに本研究の目的は高精度な流れ解析法の開発とその解析法に基づく数値解析コードの検証を行うこと、更に翼列等の実形状の流れ解析を行うことにある。
本研究報告ではまず高精度流れ解析法の根幹となる高精度化手法について述べ、数値解析コードのコアであるマトリックスソルバについて述べる。次にこの解析法に基づき作成した数値解析コードの検証計算結果について報告する。
著者らが採用した高精度化手法は、連続の式の次回残差をゼロにするために圧力のポアソン方程式のソース項の圧力勾配値を修正して新たな圧力場を求め、再度運動量保存則を解くことにより、逐次的に目標精度を達成するというものである。
数値解析コードのマトリックスソルバとして当機構が開発した残差切除法を使用している。この解法は残差方程式をＳＯＲやＡＤＩ法のごく少ない反復により求めた近似摂動量を、過去の摂動量との線形結合により最適化する方法である。に圧力場を修正することを考える。又、各方程式は無次元一般座

数値解析コードの検証は非圧縮性流れ、直交格子により平行平板間のポアズイユ流れ、平板翼列、正方形キャビティで行った。
検証結果は良好で高精度化手法の有効性を確認した。

uu    u  1   uk  p    u j  xt
x j Re  xk  xk  3 x  xk  

２．高精度化手法
流れの高精度解析を実現するための高精度化手法の論理と手順について以下に述べる。

標系表示の方程式を使用しているが、ここでは x, y, z 座標系の方程式により説明を行う。座標軸と速度をそれぞれ

x , u    1, 2, 3 とすれば x 方向の運動量保存則は（2-1）
式となる。

u
u
p     u  1   uk   u j   tx jx Re  xk  xk  3 x  xk  ・・
（2-1）
今、反復 n 回目の近似速度と次式となる。

 u その他の項
t・・（2-2）

この式の両辺の発散、すなわち div 2-1 高精度化手法の論理アソン方程式が得られる。

本手法で用いる基礎式は、運動量保存則、圧力のポワソン式および連続の式である。
運動量保存則を解いて求めた近似速度場から連続の式の残差 C  p   u     x  x t  x  xを計算し、この残差 C を 0 とする速度場 u を求めるためu  un から圧力勾配を求める

をとれば圧力のポxその他の項 ・・（2-3-a）Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 B10-1圧力勾配を求め手順 1）にもどる。

  p old   x  x ・・（2-3-b）
開始* の項は非圧縮の連続の式になっている。

Ａ○ は運動量保存則の収束に応じた誤差。

 p 
入力：  x  OL D

運動量保存則（2-1）式から速度 un p、（2-2）式から  x n が得られcal

出力： u
入力： u

連続の式る。運動量保存則が高精度に解かれていれば再計算した圧力勾配

p oldは（2-1）
式に使用した圧力勾配
x

出力：C（連続式の残差）

終了

C：と近似的に等しくなる。

OLD p   p 
入力：  x    x    

この差を  とすれば（2-3-b）が成り立つ。
（2-3-b）の圧力のポアソン方程式と運動量保存則を連立して解いても収束精度の範囲内で圧力場と速度場は変更されない。

u
t

圧力のポアソン方程式

Ａ○（2-1）式の速度 u n から計算した連続の式の残差 C は 0とならない。

 p 
出力：  x  NEW

図 2-1 時間ステップ内の計算フローチャートuC 0x

・・（2-4）
３．マ ト リ ッ ク ス ソ ル バ [残 差 切 除 法 ]（2-3）式の * で示した項は速度の発散、すなわち非圧縮性流れの連続の式である。そこで * の項を 0 となるよう強制して求めた圧力勾配値をポアソン方程式のソース項に使用する。
（2-3）式のソース項を修正した圧力のポアソン方程式は、

  P   P old    u n          ・・（2-5）
x   x   x   x    t  x   

本解法は近似解 Um にどのような修正を施せば、残差が小さくなるかという考え方で反復的に解を求めていく。
3-1

残差方程式の定式化

運動量保存則や圧力のポアソン方程式などの微分方程式を離散化して得られる差分方程式はとなる。この式から求めた p を p new として、運動量保存則

A U  f・ ・ (3-1)

を解くことにより連続の式の残差 C を遂次的に 0 とすることができる。

ここでA は 係 数 マ ト リ ッ ク ス 、U は 未 知 数 、 fはソース項である。

2-2 計算手順

ある時間ステップでの計算手順を図 2-1 のフローチャートに従って説明する。以下の検証計算では、各方程式の収束判定値ε今 、式 (3-1) の 近 似 解 を U m 、残 差 を

r m 、摂 動 とすればとしては10 13 を与えている。

量を1）圧力と速度の初期値と境界条件を与える。計算開始時の速度過r m  f  A U m・ ・ (3-2 )U  U m ・ ・ (3-3)

去値 uold

速度の初期値は連続の式を満たしているものとする。
（例えば、全領域の速度を 0 とする）
2）運動量保存則の収束解 u を後で述べるマトリックスソルバ[残差切除法] により求める。
3）連続の式の残差 C を求める。
4）連続の式の収束判定を満せば終了する。未収束の時は次の手順に進む。と 書 け 、 式 (3-1) 、 (3-2) 、 (3-3) か ら 残 差 方 程 式(3-4) が 得 ら れ る 。

A   r m・ ・ (3-4 )5）圧力のポアソン方程式を解き、圧力 p を求め、差分によりCopyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 B10-12Lr m 1     r m  1 A  m    l A   m l 1 i j kl 2

U に関する境 界 条 件 は 次 式 で 与 え ら れ る 。

U  UB： 第一種境界U  U  ： 第二種境界
l  l  B

Um・ ・ (3-5 )が境界条件を満たしていれば： 第一種境界0l3- 2： 第二種境界が 定 ま れ ば 式 (3-7) か らmU m 1 が 求 め ら れる 。残 差 方 程 式 (3-4) に 戻 っ て こ の プ ロ セ ス を 繰 り

返 す こ と に よ っ て 、残 差 の ノ ル ム を 0 又 は 最 小 に す る
解 U に収束する。

に関する境界条件は以下のようになる。

 0

さ て 、 1 ,  2   Lが 、 式 (3-8) か ら 新 し い 近 似 解

B は境界上の値であることを示す

添字

又、近似解・ ・ (3-10)








・ ・ (3-6 )残差切除法のアルゴリズム

最 初 に 残 差 方 程 式 (3-4) を ADI 法 な ど の 数 回 の反復で近似解  mを 求 め る 。次 に 残 差 ノ ル ム を 最 小とする合成摂動量

 m を (3-7)式 か ら 、 新 し い 近 似

解 U m 1 を (3-8)式 か ら 求 め る 。
L m  1 m    l m l 1・ ・ (3 -7)

l 2

図 3-1

Um 1

 U m3-3

ここでl( l  1, 2, 3, , L ) は 残 差 最 小 化 係 数 で 、l ご と の 定数である。 m は 残 差 ノ ル ムr m 1を最小とする合成摂動量である。
新しい近似解Um 1の残差r f  A Urr m  f  A U m を 求 め る 。

3） A    r m を ADI 法 等 の 最 小 単 位 の 反 復 で 近 似 解 m を求める。

m 1L r  1 A     l A  m

本解法の計算の手順は次のようになる。
（ 図 3-1 参照）
1） 初 期 近 似 値 を 与 え る 。

は 式 (3-2) 、

4） 残 差 最 小 化 係 数  l を 計 算 す る 。

 f  A  U m  A mm

計算手順2） 残 差m 1(3-7)、 (3-8) か ら 次 の よ う に 与 え ら れ る 。
m 1

残差切除法のフローチャート・・ (3-8 )ml 2Lm  l 1・ ・ (3-9 )よ っ て 残 差 ノ ル ム は 次 式 で 与 え ら れ 、 l5）  m  1 m    l m  l 1 ,l 2U m 1  U m   m を 求 め る 。

こ れ で 一 巡 し 、 2)～ 5)の 手 順 を 繰 り 返 す こ と に よ って   0 ， r m 1  0 ， U m  U  に 、又 ノ イ マ ン 問 題は最小二乗法により求める。
Copyright © 2012 by JSFM

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

では 0 ， r m 1
 は取り得る最小値に、 Um は0 乗⇒
残差ノルムU に収束する。

流れ‐15 乗⇒４. 検証結果

本数値解析コードの計算例として平行平板間のポアズイユ流れと平板翼列流れ、キャビティ流れを選び精度検証を行った。表 4-1に本検証で使用したレイノルズ数を示した。

ケース1234

表 4-1 レイノルズ数
レイノルズ数代表長さ mm3.4510.210102.6201026.9100

図 4-3 連続の式の残差履歴
横軸：時間ステップ内の反復回数

1.0

代表速度 mm/s
101550150※

密度、粘性係数は標準状態である。

※

文中ではケース 1,2,3,4 をそれぞれ Re3,Re10,Re100,Re1000 と記述している。

0.9992

図 4-4 流量の主流方向分布
青：初期 赤：収束時

ポアズイユ流れと平板翼列流れによる検証結果使用した計算格子を 図 4-1 に示した。流れはレイノルズ数 Re＝ 3 , 10 , 100 , 1000 の非圧縮性領域の流れである。

Pois DT‐4 Nm2000 NT1-8N=102E‐082×10PerrRE31.5E‐08

N=10PerrRE10Pree4-1

図 4-2 ポアズイユ流れの境界設定
青：平板壁面

1E‐08PerrRE100
PerrRE1000

5E‐09N=100-1000

012N=1500(a) 等間隔格子
101×21×21(b) 不等間隔格子
131×31×31
図 4-1 計算格子

345支流方向位置

図 4-5 圧力分布の逐次反復回数 N による変化

図 4-6 圧力の計算値と厳密解の差

4-1-1 平行平板間のポアズイユ流れPois DT‐4 Nm2000 NT100
012345‐2E‐13‐4E‐13Uerr

ポアズイユ流れの境界設定を 図 4-2 に示した。上下の青色部が平板壁面、左右の境界はすべり壁である。入口速度分布はポアズイユ流れの厳密解、出口速度はノイマン条件を与える。又、速度の過去値としてポアズイユ流れの厳密解を与え、圧力の初期値を全領域でゼロとした。よってこの計算は時間反復なしに厳密解に収束するべきものである。
図 4-3 に連続の式の残差履歴を、 図 4-4 に初回反復時と収束時の流量の主流方向分布を示した。残差ノルムは単調に減少しており約 1300 反復で収束している。
無次元流量は初回反復時には入口で 1.0、下流では約 0.999 であるが、収束時には入口から出口まで 1.0 となっている。
圧力の主流方向分布の反復回数による変化を 図 4-5 に示した。
反復が進むにつれて入口出口圧力差が大きくなる。図 4-6 には、収束した圧力計算値と厳密解の差を示した。レイノルズ数が変化しても誤差分布に差は見られない。誤差は最大 10－8 である。
図 4-7 に平板間の主流速度分布を示した。グラフは入口から出口までの速度分布を重ね書きしてある。計算した速度と厳密解の差を図 4-8 に示した。両者の差は 1.2×10-12 以下である。

主流速度‐6E‐13‐8E‐13UerrRE10
UerrRE100‐1E‐12‐1.2E‐12-121.2×10

図 4-7 収束時の平板間速度分布

UerrRE3UerrRE1000
支流方向位置

図 4-8 速度の計算値と厳密解の差

4-1-2 平板翼列の流れ

平板翼列流れの境界設定を 図 4-9 に示した。青色の部分が翼及び壁面である。
入口速度は一様流、出口はノイマン条件である。
流入部上下境界と左右の境界はすべり壁条件である。速度の過去値は一様流、圧力の初期値は 0 を与えている。タイムステップ
は 0.01 でレイノルズ数 Re 100 である。
まず、タイムステップ内の計算の収束性についてみると、図4-10 のように連続の式の残差ノルムは単調減少して約 4000 反復で収束している。初期残差ノルムは 10＋１、 収 束 時 に は 10－９Copyright © 2012 by JSFMである。
次に収束時の流量の分布図 4-11 を見てみると初期には入口から上下壁の開始位置までは 1.0、以後前縁（LE）までは 0.985、後縁（TE）までは 0.97、以後出口まで 0.985 となっている。収束時には前、後縁にわずかに流量の少ない部分がみられるものの他の断面では 1.0 となっている。
図 4-12(a) に 反 復 回 数 １ の 主 流 方 向 の 圧 力 分 布を示した。流路高さ中央の翼間の圧力を重ね書きしてある。翼の前縁、後縁にピークが表れている。ピークは翼面を含む境界面内にある。この時の入口出口 の 無 次 元 圧 力 差 は 400 で あ る 。同 図 (b) の 収 束 時の の 圧 力 分 布 を 見 る と 入 口 出 口 の 圧 力 差 は 1.0 であり、ピークも滑らかになっている。

1.0
残差ノルム-1210

図 4-9 平板翼列の境界設定
青：翼及び壁面
赤：すべり壁

図 4-10 連続の式の残差履歴

4-3

キャビティ流れによる検証

次に正方形キャビティ内の流れによる検証について述べる。計
算条件を表 4-2 に、キャビティの模式図と計算格子を図 4-14 に示した。キャビティは模式図の上面が左から右へ無次元連度 1.0で移動している。下面及び左右の壁は固定壁であり、紙面の垂直な方向に溝状の形をしている。溝の両端にはすべり壁条件を与えている｡
計算開始時の速度過去値は 0、速度初期値は上面＝1.0、その他は 0。圧力の初期値は全領域で 0 とした。
時間ステップを進めるにつれて上面に引きずられ内部でも除々に渦が形成される。
図 4-15 に連続の式の残差のノルムおよび最大値の時間履歴を示した。残差はいずれも単調に減少し、80 ステップで収束している。
文献 4 との比較を図 4-16、17 に示した。速度ベクトル、流線を見ると両者は定性的に一致している。
表 4-2 計算条件（無次元）
正方形の一辺の長さ1.0
流速1.0
レイノルズ数10, 100, 1000
格子点数
21× 21×21101×101× 5
時間刻み
0.01

上下壁の開始位置
1.0LE

第 26 回数値流体力学シンポジウム
講演番号が 時 間 反 復 １ 、(b) が 収 束 時 で あ る 。 (a) で は 翼 間の 主 流 速 度 が わ ず か に 増 速 し て い る 。(ｂ ) で は 境 界層が発達しその影響により主流速度が滑らかになっている。

TE0.96

図 4-11 流量の主流方向分布

0 ⇒TELE01.0

図 4-14 キャビティと計算格子 21×21×21400⇒(b) 収束時(a) 反復回数 1
図 4-12 圧力の主流方向分布10＋1⇒
ノルムと最大値10－１２⇒

図 4-15 連続の式の残差履歴 Re 10
横軸：時間ステップ数、Δｔ＝0.01(a) 反復回数 1(b) 収束時
図 4-13 主流速度の翼間分布

速 度 の 翼 間 分 布 を 図 4-13(a) (b) に 示 し た 。 (a)Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 B10-1
５．まとめ

高精度化手法の有効性を確認するために、平行平板間ポアズイ
ユ流れ等について非圧縮性流れ、直交格子により検証計算を行った。
検証結果から以下の点がわかった。
文献４
検証結果
図 4-16 検証結果と文献の比較 速度ベクトル図
Re 10 、格子点数 21×21×21

１．連続の式の残差の収束は単調で停留後は再び残差が増加することはなかった。
２．時間ステップを進めない定常計算においてポアズイユ流れで収束時の残差ノルムは 10－１３ となり 、平板翼列流れでは10－９となった。又、ポアズイユ流れの解は厳密解と一致した。

文献４
検証結果
図 4-17 検証結果と文献の比較 流線図
Re 10 、格子点数 21×21×21

次に断面内点数 101 × 101 、溝幅方向 5 にした格子を使って Re100、1000 の計算を行った。連続の式の残差は順調に減少している。結果の速度ベクトルと流線図を図 4-19、 20 に示した。
結果は定性的に見て不自然なところはない。

1.0E+00

Cnorm

1.0E‐020200

400

４．よって高精度化手法の有効性を確認できた。
今後は対象を非直交格子に拡張し検証を進めていく予定である。
本研究にご協力いただいたＪＡＸＡ研究開発本部の賀澤 順一
研究員、エー・エス・アイ総研の石山 毅氏の両氏に感謝いたします。

Cavity Grid101 Re=1000 Nm=2010＋0⇒

３．時間ステップを進める非定常計算により平板翼列とキャビ ty
内流れの検証を行った。平板翼列では求められた圧力、速度分布は滑らかで妥当なはがれ場にになっている。キャビティ
流れでは時間ステップ収束時の残差ノルムは 10－11 であり、レイノルズ数１０の計算結果は文献とも定性的に一致し妥当であった。

800
1000 1200
ノルムと最大値

600

1.0E‐04
norm
1.0E‐06
1.0E‐08

10－１0⇒1.0E‐10Time Step

図 4-18 連続の式の残差履歴 Re 100
横軸：時間ステップ数、Δｔ＝0.01(b) 流線図(a) 速度ベクトル図
図 4-19 検証結果 Re 100、格子点数 101×101×5●参考文献(1)田村敦宏・菊地一雄・高橋匡康, だ円形境界値問題の数値解法－残さ切除法について,日本機械学会論文集(B 編),62 巻
604 号,(1996-12), 4076-4083.
(2)Atsuhiro Tamura, Kazuo Kikuchi, and Tadayasu Takahashi.
Residual Cutting Method for Elliptic Boundary Value Problems,
Journal of Computational Physics 137,247-264(1997)Article
No.CP975807.
(3)申痢録・井小萩利明・大宮司久明, 定常非圧縮性ナビエ・ス
トークス方程式の曲線座標格子に適用できる効率的差分スキーム,日本機械学会論文集(B 編),56 巻 525 号,(1990-5),
89-1022A.
(4)川村哲也・川原睦人・平野廣和・登坂宜好・池川昌弘, 非圧
縮流体解析,東京大学出版会, 数値流体力学シリーズ
１,(1995).
(5) 菊地一雄・西澤敏雄・関根義人，残差切除法を用いた高精度流れ解析手法について，第 44 回流体力学講演会，航空宇宙数値シミュレーション技術シンポジウム 2012，講演集.
(6) 菊地一雄・西澤敏雄・関根義人，残差切除法を用いた高精度流れ解析手法について，JAXA SP 刊行予定.

(a) 速度ベクトル図(b) 流線図
図 4-20 検証結果 Re 1000、格子点数 101×101×5
Copyright © 2012 by JSFM