Li 本応用数理学会論. 文誌 Vol.8No.2.1998pp.287 305 特異な非対称行列を係数とする連立 1 次方程式での各種反復解法の収束特性と流れの数値シミュレシーョンへの応用 堀端康善 株 ) 東芝研究開発センター Solut}on of a SinglllarSystelnof Lineaτ Equations with a Nonsymmetric Coe 田 cient Matrix Using IterativeMethods and it8applica.tion to NulnericalSimula.tion of FluidFlow YasuyoshiHoribata Research& Developm 副 1t Center Toshiba Corporation Abstract.The consiste 皿 cy of a. singular system of h 皿 ea. r equations a. re discussedwhere the coe 缶 cient matrix isllollsymmetr {c a.nd the convergence propeltics of the bicon conjuga jugategradient.te gradient squa.red a 皿 d conjugate residual methods are dis cussed.numerical experiments a. re presented fbra La. pla.ce equa.tion wi もh Neumann boundaly conditions Next the consistency theory of a. sillgular system is applied to numellcal simulatioll of 且 uid 且 ow In order to ensure convelgence of the solution of the pressure equation thispaper proposes the perturbation remqva.i method using the eigenvector of the tra. nsposed coe 冊 cient ma.trixcorresponding to the eige 皿 va lueo. Therma 亅 convection of a Boussinesq Huidiss 三 mulased in a square ca vit y.the pressure equation isso 孟 ved using the proposed method.the convergellee ra.teg.of the iterative methods axe compared.the simulation result isin good a greeernent with a.be 皿 chma. rk solut 三〇 n 1 緒言 特異行列を係数とする連立 1 次方程式. tl c b の整合性と各種反復法による数値解について議論する.A が対称行列である場合について 共役勾配 法の数値解の収束特性は既に議論されている [11 ]. また 全周ノイマン型境界条件をもつボアソン方 程式を離散化してできた 特異な対称行列を係数とする連立 1 次方程式にっいて 各種反復解法で数値実験を行い 収束特性が比較されている [ 16 ]. 本論文は A が非対称行列の場合について論じる. 次に 流れの数値シミュレーションを行うときにあらわれる 圧力についての方程式との関連について考察する. 第 2 章において 特異な非対称行列を係数とする連立 1 次方程式の整合性にっいて論じる. 第 3 章 では 双対共役勾配法 自乗共役勾配法 共役残差法をこの連立 1 次方程式に適用したときの解の収 一 131 一
288 凵本応用数理学会論文誌 Vol.8No.21998 束特性を検討する. そして 第 4 章において 簡単なラプラス方程式で数値実験を行い それらの反 復解法の収束特 1 生を実証する. また 比較のため SOR 法についても数値実験する. 最後に 第 5 章 において 流れの数値シミュレーションとの関連について論じ 特異な非対称行列を係数とする連立 1 次方程式の整合化法が 従来使用されて来た非斉次項修正法よりも直接的 簡便 正確な方法であることを示す. そして この整合化法を使って 2 次元正方形キャビティ内自然対流の数値シミュレー ションし また各種反復解法の収束特性を比較する. 2 特異な連立 1 次方程式の整合性 連立 1 次方程式 2.1 ) Ax b について考える. ここで A Rnxn は特異な非対称行列である.Rnxn は実数を成分とする n n 行列の線型空間をあらわす行列 A の階数は n 1 で 固有値 0 はその特性方程式 det.xl 一 A ) 0 の単根であるとする. 但し 1 は単位行列である. 周知のように この連立 1 次方程式が解を持っための必要十分条件は 係数行列 A の階数と拡大係数行列 A A b ) とが一致することである [ 14 ]. A とその転置行列 AT の零空間を使うと この必要十分条件を次のように表現できる. A とその転置行列 AT の零空間をそれぞれ N A ) N AT ) とすると 1>A ) ; span { e } N AT ) ; pan { * } と書ける. ここで e と は それぞれ A と AT の固有値 0 に対応する固有ベクトルである す ると A の値域 R A ) は N A : ) の直交補空間になる. すなわち 実数を成分とする n 項列ベクト ルの線型空間 Rn は部分空間 R A ) と N A τ ) の直和である. Rn R A ) 千 N AT ) そして 連立 1 次方程式は b のときにのみ解を持つ 同様に 随伴方程式 AT π b はゲ e のときにのみ解を持つ. 仮定により固有値 0 が待院方程式 det a;1 A ) 0 の単根であるから 行列 A のジョルダン標準 形は 固有値 0 に対するジョルダン細胞として O 1 ) 0) を一っだけもつ. したがって 行列 AA の階数は行列 A の階数と等しくなり となる. 7 α nk : AA )? αnk A ) n 1 固有ベクトル e と について次の関係が得られる. 定理 [1] 固有ベクトル e と e * は直交しない. e e * ) 0 一 132 一
特異な非対称行列を係数とする連立 1 次方程式での各種反復解法の収束特性 289 証明 : e e ) 0 を仮定す る. すると A の値域 R A ) は N AT ) の直交補空間であるから e R A ) となる.N A ) の基底く e > を拡大して Rn の基底く e e1 e2.en _ 1 > を得たとする. このとき < AeiAe2 _ Aen _ i > が R A ) の基底になるから e は Ae1 Ae2.Aen_ 1 の線型結合として表わされる これより 崩 c el + c2e2 + + ctt 一ユ r. ) A e c1 /1e1 + c2ae2 + _ + c _ iae L 1 A + c2ae2 +...+ c. Ae. ) A となる. 一方 んAe o であるから 互いに線型独立な Cle ユ十 C2e2 +... 一 Cn _ 1en _ 1 と e は行列 AA の零空間のベクトルになる. したがって 行列 AA の値域 R AA ) の次元は n 2 以下となり 結局 r α 淵舶 ) n 2 となる. これは矛盾である. 口さらに 次の定理を得る. 定理 [2 ] 線型空間 Rn は部分空間 N A ) と R A ) の直和である. Rn N A ) 千 R.A ) 証明 : R.A ) の基底をく i f2... i _ 1 > とする.e fliず 2...fn_ 1 の間に線型関係 C e + Clfl + C2 /2 +.. + C 一 fn 1 があれば これと e * との内積をとると CQ e e ) 0 を得る. 定理 [1 コより e e * ) 0 であるか ら Co 0 となる. したがって 線型関係は Clf + C2 / 2 +. + Cn fn 1 となる. < f1 2_ fn_ 1 > は R A ) の基底だから Cl c2. Cn 1 0 となる. したがっ て ベクトル e f1 f2.. fn 1 は線型独立になり Rn の基底になる. ロ 3 各種反復法の収束特性 反復法を使って特異な連立 1 次方程式 2. 1 ) を解く場合の収束特性にっいて検討する. 3.1 双対共役勾配法 BCG 法 ) 双対共役勾配法を特異な連立 1 次方程式 2.1 ) に適用する. 双対共役勾配法では ムにしたがってベクトル列 Xl x2 _ を求める [6 ]. o と x6 は初期値である. 次のアルゴリズ アルゴリズ厶 1 ro b IXo l r6 b 一 ATxe Po ro : P δ rg ; k : O while IFrkil Odo α k rk rx.)/ APk P :. ) 鰥 +1 Mk + 叫 P た rk rk 一 drkapk ; 噺 1 β た rle +1 rz + 1)1 rk r 蒸 ) 二 rx. 一 katpr Pk + 1 rk + 1 + fikpkl px. + 1 丁誕 + 1 + βkp 壽義 :kt 斗.1 一 133 一
290 日本応用数理学会論文誌 Vol.8Na21998 このアルゴリズムによって得られるベクトル列 {Pi } と { ri } に対して双対直交性 と共役直交性 丁凱丁ゴ ) o i ゴ ) が成り立っ [ 61. P 菷 4P ゴ ) 0 i ; ゴ ) 定理匚 2 ] から 連立 1 次方程式の解 x は次のように書ける. x y 十 z 但し y 1> A ) z R A ) ここで z は一意的に決まるが y は一意的に決まらない. b R A ) ならば ri R A ) となる. そして 探索ベクトル Pi も R A ) のベクトルであるから h 1 次元の部分空間 R A ) の中で解を探索し z を求めることになる. 定理 [2 ] より e と Pt は線型 独立であるから ベクトル e は探索過程に参入しない. 探索が続く限り 双対直交性から ro rl... ri は線型独立になる. そして 定理 [2 ] より e と ri は線型独立であるから e r rl.. 物は線型独立になる. したがって ある m n 1 で rm o となって解が求まる. 32 自乗共役勾配法 CGS 法 ) 自乗共役勾配法は双対共役勾配法におい 1 r6 re としたものである. 次のアルゴリズムにした がってベクトル列 M1 M2.. を求める [ 15 ]. o は初期値である. アルゴリズム 2 ro b Amo PO ro i k O while Ilrkil @O@ α 沁 ; ア和 痢 / APk r j hk + 1 ε た一 ak Ap ) XA + 1 需鞠 + 蝋 e 杜幅 ) 併婦二恥一蝋 e ん嗣 ) β k ro ra: + 1 )1 r Crk ) ex + 1 rk + 1 + 距んみ +1 Pk + 1ek+1 + β A hk+1 βkpk ) 斥 ; k : 十 1 共役勾配法の場合と同様に b R A ) ならば ある m n 1 で甲 η o と って解が求 NII-Electronic る. 一 134 一 N 工工一 Eleotronio Library
特異な非対称行列を係数とする連立 1 次方程式での各種反復解法の収束特性 291 3.3 共役残差法 CR 法 ) 共役残差法では 次のアルゴリズムにしたがってベクトル列 Xl x2 _ を求める [ 13 ].Xo は初期値である. アルゴリズム 3 ro b Axo Po ro ; KlO while l1 副 Odo α k rk APk ) 1 Apk APk ) 記軒 1 コ Ck + α kpic rk + 1 rk akapk βk 一 Ark + 1 APk )/ APk Apk ) Pk + 1 Tk + 1 + βk Pk k k 十 1 b R A ) ならば ri Pi R A ) となる. したがって 部分空 P7 R A ) の中で解を探索することになり ベクトル e は探索過程に参入しない. この場合 rk APk ) ric Artl ) 0 である限り 探 索を続けることができ rb ri ) の値は単調に減少し続ける. そして 解の成分 z R A ) が求めら れる. 4 ラプラス方程式での数値実験 4.1 連立 1 次方程式 次のラプラス方程式を考えることにする. 解領域は Ω 0 1 ) 2 である. 2 一 i11 Ω 象一 n Ω ここで g は u x ) 二 Ml + x2 が解になるように決める. 格子点が境界付近に集中するように次式に したがって座標囮 x2 ) を座標 ξ1 ξ2 ) へ変換する. β+ 1 ) [β + 1 )/β 一 1 2 ) 亅 ξ1 一 β + 1 1 2 { 1 1 Kβ 一ト 1 )/β 一. 1 ) 2 ξ1 1 } β + 1 )[ β + 1 )/β 一 2 ].) 亅 ξ? 一 β+ 1. 2 : 2 { 一 β 一ト 1 ) ノ β 一 1 ) ] 2 一 ξ2 1 } ここで β はパラメータであり β 認 Ll にとることにする. この座標変換によって 解領域 Ω を写像してできた計算領域 Ω O 1 ) 2 上で離散化する. 格子間隔は 1 / 29 にとる. 特異な非対称行列 A R9 OOxgoo を係数とする連立 1 次方程式. 勉 b が得られる. しかしながら この連立 1 次方程 式は擾乱を含み b R A ) となる. 得られた連立 1 次方程式から 次のようにして擾乱を取り除く. 1. 行列 A の転置行列の零空間 N.4T ) のベクトル e を方程式 ATe * o から求める. ベクトル の一つの成分の値を固定してから反復解法を適用すれば容易に求まる. 一 135 一
292 日木応用数理学会論文誌 Vol 8No.21998 2. 行列 A の値域 R A ) は N AT ) span { e * } の直交補空間であるから 方程式 Am のベクトル b を R A ) へ射影すると b の右辺 が得られる.b をこの br で置き換える. b 屠卜 b e * )IIIe lj2 ] e * こういうふうにして擾乱を取り除いた連立 1 次方程式 Ax br は br R A ) となり 解をもっ. 擾乱の影響を考察するときは この連立 1 方程式に擾乱 b を加える. そのときの方程式は次の になる. よう Az b ここで b br 十 b 擾乱 b は b r!lle il) e の形で与える. ここで ty の値は 二 7 ]lbrll δ 1 > / i : 歪 0 く δく 1 ) から決める.b と R A ) の間の角度を θ とすると この δ について の関係がある sin θ ilbll11ibll δ 42 数 { 直実験 CRAY Y MP C90 1GFLOPS ) を使う. 擾乱が反復法の収束特性に及ぼす影響を広範囲かつ詳細に考察するため 倍精度実数型の浮動小数点数の演算 計算機イプシロンは約 2.5244 lo } 29 ) で行う. 双対共役勾配法 自乗共役勾配法 共役残差法 そして比較のため SQR 法 一 136 一
特異な非対称行列を係数とする me i fl1 次方程式での各種反復解法の収束特性 293 4 コ 丶 ρ 1. リ < ). b 一 0 一 4 一 8 一 12 一 16.20 一 24 0 15 30 45 number of iterations 60 75 Fig.1.PBCG fordifferent perturbatio 皿 s. を加えた 全部で 4 種類の反復解法で擾乱のない場合と擾乱を含んだ場合を解く. 初期値は常に Xo. 褐 o にとる. 但し 双対共役勾配法自乗共役勾配法 共役残差法にっいては 収束を加速するため前処理を施す. 前処理には不完全 LDU 分解を使う. すなわち 行列 A を不完全 LDU 分解し て LDU の形にする. 左下三角行列 L と右上三角行列 σ の非対角要素は元の行列 A の対応する要 素に等しくとり 対角行列 D の対角要素だけを計算する [13 亅. そして方程式 LD 一 1 の.4 ゐD 冖 lb の に双対共役勾配法 自乗共役勾配法 共役残差法を適用する. 前処理っき双対共役勾配法 PBCG 法 ) で方程式を解いたときの結果を Fig.1 に示す. 擾乱を除い た方程式では Xi は解に収束する 擾乱があると Xi. はある点まで収束するが途中から残差が増大 し始め っいには発散してしまう. そして 擾乱の大きさが大きいほど早い段階から残差が増大し 始める. これは 次のように説明される. 擾乱が存在するため rk r 毒 ) はある程度以下には小さく ならない. 一方 探索ベクトル Pic と域にそれぞれ擾乱として e とげの成分が入っていようと eg /> A ) N AT ) であるから 五 p 齢露 ) の方は rk r 壽 ) に構わずに小さくなり続ける. その ため ある地点から ak rk r X )1APk pi ) の大きさが擾乱の無い場合よりも大きくなり始め 発散し始める. そして 擾乱の大きさが大きいほど rk の 噌最小値は大きくなるから より早い段階 から発散し始める. 前処理っき自乗共役勾配法 PCGS 法 ) で方程式を解いたときの結果を F 三 g.2 に示す. 前処理っき双対共役勾配法と同様の傾向が得られているが 前処理つき自乗共役勾配法の方が収束速度が速い. また 同じ大きさの擾乱では前処理つき自乗共役勾配法の方がより小さい残差になってから発散を始 める. 前処理っき共役残差法 PCR 法 ) で方程式を解いたときの結果を Fig.3 に示す. 擾舌 L が含まれてい ない場合には Xi は解に収束する. 擾乱が混入している場合には mi は擾乱のない場合と同様にある 点まで収束し続ける. しかし ある点からは反復を続けても残差の大きさはほとんど変化しなくなる. そして 擾乱の大きさが大きいほど より少ない反復回数でより大きい残差に飽和する. これは次の ように説明できる. mi が方程式 Ax br の解に近づくと 残差はハ b 4 鞠 b 酬 le 類 ) に すなわち行列 AT の零空間のベクトルに近づく. そのため rxt Apk ) ATrk ph ) の値は 0 に 近づき )k の値は小さくなり 反復解の改良は止まってしまう. 擾乱の大きさが大きいほど より早い段階で 残差の中で行列 AT の零空間のベクトルが支配的になり 探索が実質的に止まってしまう. 一 137. 一
294 )] 本応用数理学会論文誌 Vol.8 No.21998 4 一 丶 ρ < ) o[ bdo 0 一 4 一 8 一 12 一 16 一 20 一 24 0 0 20 30 40 number of iterations 50 60 Fig.2 PCGS fordifferent perturbati 皿 s. 0 一 5 丶 皿 の ーガ < ) o[ bdo 一 一 5 一 10 一 15 一 20 一 25 0 40 80 120 160 200 number of iterations 240 280 Fig3.PGR.fordifferentperturba. tions. SOR 法で方程式を解いたときの結果を Flg.4 に示す. 加速係数 ω の値としては最適値 ω pt 190 を使用した. 前処理つき共役残差法と同様の傾向が得られている. 次に 擾乱のない方程式について係数行列の特異性を消してから反復解法を適用する. 特異性を消すために 解 x の第 2 成分 :n2 の値を固定して :; 2 O とし 連立 1 次方程式 の第 2 式をこれと入 れ替える そういうふうにして得られた連立 1 次方程式は A x bt となり 係数行列 A は正則である. この連立 1 次方程式を前処理つき双対共役勾配法 前処理っき自乗共役勾配法 前処理つき共役残差法 SOR 法 ω pt 1.98 ) で解き 元の特異な連立 1 次方程式を そのまま解いた場合と比較する. この場合も初期値は Xe ; 端 o にとる. 結果をそれぞれ Fig.5 Fig.6Fig.7Fig.8 に示す. いずれの場合も 特異な連立 1 次方程式をそのまま解く方が 特異性を消去してから解くよりも収束が速い. 一 138 一
. 特異な非対称行列を係数とする連立 1 次方程式での各種反復解法の収朿特性 295 20246802468024 の 丶 ρ < ) o 冖 boo [ 一幣一一一璽 葺 111222 層一一一 0.00 200 300 400 500 number of iterations 600 Fig4.SOR wopt 1.90 ) fbrdsfferent perturba.tions. 至 萋. 乱メ ) ).. b [ 三莖亀島嫋毳 8 2 0 一 2 一 4. ロ6 一 8 一 OI 0 10 20 30 40 50 number of iterations 60 Fig 5PBCG fbrthe si 皿 gula. r and nonsingular s ] stenlg. ρ ミ ρ ) 一 [ 蚕醗易 一萋 炉 尭く ) ぎ [ 薫島楞譽 8 6 4 2 0 2 4 6 8 10 0 10 20 30 40 number of iterations 50 Fig.6PCGS forthe singula. r a. nd nonsingular g. } ysterng.. 139 一
296 目本応用数理学会論文誌 Vol 8 No.21998 ρ ρ < ) 三 碧 磊 メ ). 二糞誌易 8 0 一 2 一 4 一 6 一 8 一 10 0 50 100 number of i 重 erations 150 Fig.7.PCR for the singular and 皿 onsingu ar systems. ニ 互 < ) 葛毎 蚕一 暫 o. メ ). b 莖乱 甥 8 1 0 1 2 3 4 5 6 7 8 9.10 0 50 100 150 200 mmber of iterations 250 Fig.8.SOR forthe singular and nonsingular systems. 5 流れの数値シミュレーションへの応用 5.1 特異な連立 1 次方程式の整合性と圧力についての方程式 非圧縮性流れや低マッハ数流れを数値シミュレーションするとき ノイマン型境界条件をもつ圧力 p についての方程式 一一 S in v 嘉一 G nc を解く必要がある国 このソース関数 3 と墳界値関 ikg との間にガウスの発散定理により 5 1 ) Sd v ff f の関係が成り立っ. ここで レは解の領域 σ は解の領域の境界をなす閉曲線 p/ n は C の法 線方向に沿った p の微分係数である. しかしながら 方程式を離散化した場合 離散化誤差のために Gds 一一 140 一
The Japan Society Soolety for Industrial 工 ndustrlal and Applled Applied Mathematlos Mathematics 特異な非対称行列を係数とする連立 1 次方程式での各種反復解法の収束特性 297 一般には式 5.1 ) の積分拘束条件は満たされない. そのため 離散化して得られた連立 1 次方程式を 反復解法で解くと発散してしまうことがある 2 3 7 ]. また 流れ場の中のある一点での圧力の値を 規定してから数値解を求めると発散は抑えられるが しばしば解が歪んでしまう [9 10 ]. そのため 通常 非斉次項修正法が使われる [1 2 3 7 91. 非斉次項修正法は ソース 関 iks を 次のように修正する. s * s / v ここで ε 一 ffv s 丿 fc dv G ds S を S * で置き換えると 積分拘束条件は正確に満たされる. しかしながら E を求めるための面積分と線積分を精度よく数値積分するのは非常に煩雑である [ 2 9 ] この問題を特異な連立 1 次方程式の整合陵という観点からみると次のようになる. ノイマン型境界 条件をもっ圧力についての方程式を離散化すると 特異行列を係数とする連立 1 次方程式になる. 式 5.1 ) の積分拘束条件はこの連立 1 次方程式の整合性条件と等価である. この連立 1 次方程式を反復解 法で解いたときに発散するのは 前章で考察したように 連立 1 次方程式の中に擾乱が混入している ためである. また 流れ場の中の基準点での圧力の値を固定してから数値解を求めると 流れ場に基 準点の位置についての依存性がみられるのも 連立 1 次方程式の中に混入している擾乱のためである. 物理空間で等間隔の格子点を使って離散化すると係数行列は対称になるが 一般座標へ座標変換し てから計算空間で等間隔の格子点を使っ て離散化すると係数行列は通常非対称になる. 非対称の場合 には 前章までにおいて論じたように 特異な非対称行列を係数にもつ連立正次方程式を その係数 行列の転置行列の 固有値 0 に対応する固有ベクトルを使って擾乱を取り除いて整合化させると発散 が抑えられる. また 擾乱を取り除くと 基準点で圧力の値を固定しても基準点依存性はなくなる. この連立 1 次方程式の整合化法の方が非斉次項修正法よりも直接的であり である. そしてより簡便で 正確 5.2 2 次元正方形キャビティ内自然対流での数値実験 5.2 1 支配方程式 Fig9 に示した 2 次元正方形キャビティ内自然対流について考える. 上下の壁は断熱壁であり 左一定で TH 305 K Tc 295 K とする. 流体のプラントル数 P7 u /X 右の壁の温度はそれぞれ を O.71 にとる. レイリー数 Ra β9α 甘一 Tc ) D3A u は 106 にとる. ここで β は 300K での体積膨張率 X と v はそれぞれ温度伝導度と動粕性係数 g は重力加速度 D はキャビティの辺の長さである. この自然対流を記述するのに 低マッハ数流れの方程式 [8 を使う. 低マッハ数流れの方程式は次 のようになる. ここでテンソル表示を使い 添字 1 ゴは.7; x 軸のベクトルおよびテンソルの成分 に対応する値 1 2 をとることにする. まず ナビエ ストークス方程式は P 霊占 Pd σ 1ぽ一 ρ 万 τ + 蕊爾 D?t Xid ; σ 2で ρ 万 τ 一 + 冨 ρ 9 一 141 一 NII-Electronic N 工工一 Eleotronlo Llbrary Library
298 凵本応用数理学会論文誌 Vol.8No.21998 Z T Z 一〇 縛 ldi 1 TH T 一 \ T Z 一〇 〆 X C ー 9 D 一一一一 > Fig.9.Thermal convection in a square cavity と書ける. ここで u w はそれぞれ速度の : z 成分 ρ は密度 Pd は圧力の動力学成分 転ゴは粘性 ストレステンソルで速度から鮪 一 諏 謙一 1 聯 とあらわされる. μ は粕陸係数である. 連続方程式は 1D6 1 Pr 砺 EDt Pr dt 賜 となる. ここで E は単位質量あたりの内部エネル ) e 一 p は流れ場の代表的な圧力である. エネル ギー保存の式は ケ 1 一離 響 ここで 7 は定圧比熱 Cp と定積比熱 Cv の比で 7 Cp / Cv である. qi は熱流束密度で温度丁と熱 伝導度 κ とから となる. 状態方程式は次のようになる. qi 一 κ T 爾 P. r ツー 1 )ρ E 圧力についての方程式はナビエ ストークス方程式から次のように得られる. 5.2 ) 1 ) 一一畠 ここで M. 塞 窪 ) 黔 黔 M 宣ひはそれぞれ Navier Stokes 方程式の z ご成分の移流項 粘性項 重力項の合計をあ らわす 境界では u ; t 疋 0 となることから ナビ エストークス方程式よりノイマン型境界条件伽 G π が得られる. 一 142 一
The Japan Society Soolety for Industrial 工 ndustrlal and Applled Applied Mathematlos Mathematics 特異な非対称行列を係数とする連立 1 次方程式での各種反復解法の収束特性 299 Fig10.Grid system. 5.2.2 数値シミュレーション 特異な非対称行列を係数にもつ連立 1 次方程式の整合化法を使って圧力にっいての方程式を解き キャビティ内の自然対流を数値シミュレー. 一ションする. 第 4 章では 擾乱が各種反復法の収束特性 に及ぼす影響を広範囲かっ詳細に考察するため rt b Axi の大きさと右辺のベクトル b の大きさの比が 10 20 倍精度実数型の浮動小数乾数の演算を採用し 残差 以下になるまで反復を実行し た δ 1]bll /Hbll 10 22 となるような擾乱 b まで調べた. しかし 流れの数値シミュレーションで は 残差 ri b AXi. の大きさと右辺のベクトル b の大きさの比が 10 1e ま 程度まで収束すれば実用 上十分である. それで CRAY Y MP C90 を使い 単精度実数型の浮動小数点数の演算 計算機イプ シロンは約 7.1054 10 15 ) で行う. 支配方程式を差分化するために 50 50 の格子点を使う. 解の精度を向上さすため 座標変換を使い Fig10に示すように格子点を境界付近に集める. 差分化により 式 5.2 ) から圧力にっいての連立 1 次方程式 Ax b が得られ る. 係数行列 A R2500 25GO は特異な非対称行列になる. 差分化した方程式を MAC 法に類似の方法 [8 亅を使い tx / D2 z 1.08 ]0 4 なる時間刻み幅 t で時間方向に積分する. 静止状態に瞬時に TH Tc の温度差を加える衝撃出発とする. 各時間ステップにおいて圧力についての連立 1 次方程式を解く必要がある. 初期の時間ステップで は この連立 1 次方程式に δ lbll /llbll re 10 5 程度の擾乱 b が混入していた. 第 4 章からわかる ように 擾乱が混入したまま反復解法を使うと 収束する前に解が発散したり 停滞したりする. そ こで 各時問ステップにおいて 整合化法を使ってこの擾乱を除去してから 双対共役勾配法 自乗共役勾配法 共役残差法 2 色 SOR 法 の 4 種類の反復解法で解く. 双対共役勾配法 自乗共役勾配法 共役残差法の前処理には 不完全 LDU 分解に Gustafsson 流の補正を行った MJLDU 分解を用いる [ 131. MJLDU 分解に用いる補正係数 α の値を O.98 とする. 各反復解法において 初期値 Xo を決めるのに 一 143 一 NII-Electronic N 工工一 Eleotronlo Llbrary Library
300 日本応用数理学会論文誌 Vol.8No.21998 400 切唱 O 喟一 邸 〇一順 隔 O O 自 5 口 003 002 001 0SOR R & B) PBCG PCR PCGS Fig.11.Average number of itera. も ions 各時間ステップにおいて 初期値 XO を mo 0 にとる 各時問ステップにおいて 一つ前の時間ステップでの圧力を初期値鞠にする の二通りの方法を採用した. そして Il Ami bii /11 b ll く 10 10 に収束するまで反復を続ける. 時間ステップ数が 10 以上になると 擾乱 b の大きさは常に δ IIbll /1 re 3 10 6 であった. どの 反復解法を使っても また初期値鞠の決め方にかかわらず 流体変数は同じ時間的経過を辿り そして時間ステップ数が約 1500 位になると定常状態に到達した. 定常状態に到達するまでの 各反復解法の 1 時間ステップ当たりの平均反復回数と平均 CPU 時間 をそれぞれ Fig.1L Fig.12 に示す. 初期値を常に Xe 数は同数になり Gllstafsso コi 流の補正を行った前処理つきの双対共役勾配法 自乗共役勾配法 共役残差法 および 2 色 SOR 法 u ) 二 ept 1.93 ) でそれぞれ 45 29 69 353 回であった. 当然各時間ス テップにおいて所要の CPU 時間はほぼ同じ長さになり Gusta.fsson o にすると 各時間ステップにおいて反復回 流の補正を行った前処理つきの 双対共役勾配法 自乗共役勾配法 共役残差法 および 2 色 SOR 法 ω pt 1.93 ) においてそれぞれ 50 3244 97ms であった. 一方一つ 前の時間ステップでの圧力を初期値にすると 時間ステッ プが進むにしたがいより少ない反復回数で収束するようになる. これは 定常状態に近づくにしたが い 一っ前の時間ステップでの圧力が連立 1 次方程式の解のより良い近似値になっていくからである. Gustafsson 流の補正を行った前処理つきの双対共役勾配法 自乗共役勾配法 共役残差法 および 2 色 SOR 法 wopt 1.93 ) での平均反復回数はそれぞれ 251326 157 回であった. また 平均 CPU 時間はそれぞれ 32 16 19 52ms であった. 初期値記 } の決め方にかかわらず 自乗共役勾配法の 所要 CPU 時間が最短になる. そして 共役残差法 双対共役勾配法 2 色 SOR 法の順番である. 5 2.3 係数行列を正則化した場合の収束特性初期値を常に Xo 蕭嬬 o にとった場合の各反復解法の収束特性を調べる. 前述したように この場合の収束特性はどの時間ステップにおいてもほぼ同一になる. したがって 定常状態に到達した後の時間ステップでの すなわち定常状態の速度分布と温度分布から決まる 圧力にっいての連立 ] 一 144 一
特異な非対称行列を係数とする連立 1 次方程式での各種反復解法の収束特性 301 100 80 傍邑 60 り葺 P 40 u 20 0SOR R & B ) PBCG PCR PCGS Fig.12.AverageCPU time 次方程式 Ax b について調べることにする. この連立 1 次方程式に混入している擾乱 8 の大きさは である. δ ; llb 闘 /1 bll;2.475 x lo 6 まず ベクトル e * を使ってこの擾乱を除去した後 Gustafsson 流の補正を行った前処理つきの双 対共役勾配法 自乗共役勾配法 共役残差法 および 2 色 SOR 法 w pt 二 1.93 ) で解く. 次に 擾乱を除いた方程式にっいては基準点での圧力を固定しても解は歪まないから 基準点で圧 力を固定して係数行列の特異性を打ち消してから反復解法を使ってみる. 第 1251 番目の成分 1251 が 定義されている格子点を基準点とし その値を 0 に固定して係数行列を正則化してから連立 1 次方程 式 A x bt を解く.Gustafsson 流の補正を行った前処理つきの双対共役勾配法 自乗共役勾配法 共役残差法 および 2 色 SOR?k w pt 1.9.3 ) ついて特異なまま扱った場合と正則化した場合の収束特陸の比較をそれぞれ Fig.13 Fig.14Fig15 Fig16に示す 第 4 章の結果と同様に いずれの場合も特異なまま扱った場合の方が収束速度が速いのがわかる. 5.2.4 数値シミュレーション結果の検証 どの反復解法を使っても また初期値をどちらの方法で決めてもほとんど同一の定常状態が得られ た. 速度分布と温度分布をそれぞれ Fig.17Fig18 に示す. 温度は? 二 [T と無次元化してある. 得られた定常状態から次のような特徴的な量を計算する. NUONUn 鉛直壁侮 0 ) 上での平均ヌッセルト数 萄 max max axnumin 鉛直壁幕侮 0 ) 上での最大局所ヌッセル ト数 及びその位置 ) 鉛直壁師 0 ) 上での最小局所ヌッセル ト数 及びその位置 ) 鉛直 中線侮 α5 ) 上での最大水平速度成分 及びその位置 ) 水平中線 写 0.5 ) 上での最大鉛直速度成分 及びその位置 ) ここで 座標と速度は次のように無次元化してある. 7 }7 + 距 )12 レ TH Tc ) 記罵 m /D 署 z /D 一 145 一
The JapanSociety Society forindustrial Industrial and andappliedmathematics Mathematics 302 H lsc fief M twlgle}ijs}:ra JSa ± Vol. 8No21998 -- - -3 {SmoN-- -4 o-oo-oooo-- ek L'si' ge es a E 1 o -1-2 -5-6 -7.8-9-10 o Fig. 13. PBCG 10 20 30 40 50 number of iterations forthe singular and nonsingll}ar systems for pressure. AF ps)dd 2 o ' rv 1.-.Xts-H-2 vy-ov e.-oo"uaoo-- -4 :L' -6 stse 1uy :- 'z -8 2-10 O 10 20 30 40 nurnber of iterations Fig. 14. PCGS for the singulat and llonsingular systems for pressure t'-s.-x D ))b9.. -2 " - ti$-voe -4 :u diss 1u 'Z 8 o -6-8 ge ff. -10 o 20 40 60 80 IOO 120 140 number of iterations 160 Fig. 15. PCR forthe singular and nonsingular systems for pressure. A-iig$i9.- e"tu-uaoo'- -146- NII-Electronic Library
特異な非対称行列を係数とする連立 1 次方程式での各種反復解法の収束特性 303 一至 ) 鳥巷 鹽 三ー メ 一 ).. 莖蕊 蒭 5 1 0 1 2 3 4 5 6 7 8 9 10 亀 點 鹽 L 9.. 曽驢 馬 辱.. 冒.r 辱巳辱 non singular 幽. singular ω P 193 ) 0 60 120 180. 240 number of iterati.ons ω P 嵩 L97 ) 300 360 Fig.16Red BlackSOR for the singulax a 皿 d nonsingular syg.tems for pressure ttt 1 鹽 一一一 } 一一一辱.. 鹽 畠 L r 一鱒 r 隔亀噛.r 幽. 巳一.. 曾 凾一一. 臨鹽 げ ノー丶 r 9 曹辱一辱嘘 凾一尸 M M 馳 膨篠 ここごここ : 二 :: 二三三三三 i ミ 1 : 翩 lll ::: :1 誉 lllll 三 i 三三淵曲 r 玩 _..._ 9. 二 :: 隔 W ::llllll 常 1::: :::::::::1 ::IIIi ttl LLIIil 曲 11..1.1..1. 鹽幽... lltll.. ili r lh..11 ; lilil. 層... 層... 脚 ::: :: 11 Ul 三藍瀧 1 譱三謡 i 鬢 ; 鱗 11h _ 一一一一一一一一一. tlttilv iill 1 三 :::::: こここごこここ三諺嬲 ; 1 こt:.:: :::::::: ; 111 : ここ 彰二 ll::. 匿隔 鼬 ttf Fig.17.Vel city vect rs. Fig.18.Isotherms.The col 且 tours ra. nge froln 0 5 to O.5 wit 1 a contour interyalof O.1. 一 147 一
304 日本応用数理学会論文誌 Vol 8No.21998 Tablel.Differences % )betweei1も he solutio 皿 reported i 皿 this paper a 皿 d a benchmark solution Nu NUmat[ @ 乏 NUmi 質 a. ax @2 O n[ 毓 Presentr 器 uits 8.706 17.799 0.0344 0.952 1 67.75 0.847 224.06 0.0374 de ValllDavis 8.8 董 7 17925 0.0378 0.989 1 64.63 0.850 219.36 G.0379 Difference %) 1.3 0.7 9 0 3 7 0 4.8 0.4 2.1 13 a up /x1 wd!x 計算値をベンチマ ーク解 [ 51 と比較すると Table1 のようになる. 誤差が最大になるのは最大局所ヌッ セルト数の位置についてであり 9.0% である.2 番目に大きいのは鉛直中線上での最大水平速度成分にっいてで一 4.8% になる.2 次元正方形キャビティ内自然対流については 既に数多くの数値シミーュレションが行われ ベンチ マーク解と比較されている [ 4 12 ]. それらの中で本論文の結果の誤差は小さい方であり したがって妥当な定常状態がシミュレーションできていると考えられる. 6 結言 旨 特異な非対称行列を係数とする連立 1 次方程式の整合性にっいて論じ 双対共役勾配法 自乗共役 勾配法 共役残差法を適用したときの収束特陸を考察した. 擾乱がなく 整合がとれているとどの反 復解法でも解が得られるこ とを示した. このことを簡単な数値実験で実証した. 擾乱を混入さした数 値実験では 解が発散したり 停滞するのが観測され その理由にっいて考察した. また 擾乱のな い連立 1 次方程式では 特異のままで解く方が 解の一成分の値を固定して係数行列を正則化してから解く場合より収束が速いのがわかった. 次に 特異な非対称行列を係数とする連立 1 次方程式の整合性と 流れの数値シミュレーションと の関連を論じ ノイマン型境界条件をもつ圧力についての方程式を安定に解くために 従来よく使わ れてきた非斉次項修正法の代わりに 係数行列の転置行列の固有値 0 に対応する固有ベクトルを使っ て連立 1 次方程式を整合化させる方法を提案した. 整合化法の方がより直接的で 簡便で 正確であ る. 整合化法を使って 2 次元正方形キャビティ内自然対流を数値シミュレーションし 得られた結果 をベンチマーク解と比較した. 良好な一致が得られ 整合化法の有効性が確認できた. そして 双対 共役勾配法 自乗共役勾配法 共役残差法 および 2 色 SOR 法の所要反復回数と所要 GPU 時間を比 較した. また 圧力についての連立 1 次方程式を特異なまま扱った方が 基準点で圧力の値を固定し て係数行列の特異性を打ち消してから反復解法を使う場合よりも収束速度が速くなるのがわかった. 謝辞 本研究中 山形大学工学部機械システム工学科賈為博士 九州大学大学院数理学研究科田端正 久教授 株 ) アールフロー竹田宏博士から有益な御教示をたまわった. 記して深謝する. 参考文献 国 Anderson D.A.Tannehil1.C. and Pletcher ComputationalFhid Mecb 細 cs and Heat Transfer Hemisphere New York1984. [2 亅 BiringenSand Cook C. On pressurebo 匸 mdary conditions for the l 皿 compressible Navie ト Stokes equations using nollstaggered ハ厂 umer.heat 7b ansfer grids 13 1988 ) 241 252. 一 148 一
特異な非対称行列を係数とする連立 1 次方程式での各種反復解法の収束特性 305 [3 亅 Briley W R.Nulnerical Inethod for predictsng th reedimensional steady viscous fiowinducts J.Comput.Phys. 14 1974 8 28. ) r41 de Vahl Davis G.and Jones 1.P Natu 凪 lconvection of aiτ in a aquare cavity : acomparison exercise nt.jl IVumer. ル1 etho オ 3 Fluitls 3 1983 ) 227 248 [ 5 ] de Vahl Davis G. Natural co 皿 vection of air i 皿 a. aquare cavity : abench mark soeution Int. ハ ume? Methods Fluids3 1983 ) 249 264 Co 皿 jugate bcrrad ient methods fbr inde 且皿 ite systems LectureNotes 506 1976 ) 73 89Springer V611ag Berlin. 16] Fletcher R Math. [ 7 ] Ghia U Ra 皿 amurti R. and Ghia K.N. Solutionof the Neumd nn pressureproblem j 皿 genetal ortbogonal coordinates using the mul 七 igridtechnique AJAA J. 26 1988 ) 538 547 8 ] Horibata Y.Numerical simulation of a low Mach number fiowwith a la e temperature va.riation Fluids 21 1992 )1185 200. [9 ] 賈為 中村佳朗 保原充 u v p 法による自然対流の数値計算法 ながれ 9 1990 ) 34 52. [10] 賈為 有限体積 差分法による非圧縮性流の計算法に関する研究 名古屋大学博士論文 1991. [ 11 ] 1aasschieter E.F. Preconditioned conjugate gradientsf{)r solving singular systems.o / Oomput.Appl.Math.24 1988 ) 265 275. [12 ] Markatos N.C.and Perjcleous K A. Laminal and turbulent natural couvection ina enclosed cavity nt. 訊 HeαtMass 7} ansfer127 1984 )755 772. [13 ] 村田健郎 名取亮 唐木幸比古 大型数値シミュレーション 第 6 章 岩波書店 東京 1990. [ 141 齋藤正彦 線型代数入門 東大出版会 東京 1969. {IS 亅 Son 皿 eveld P. CGS : a.fastla.nczos type solver {fornonsymmetric linearsystems SIAM J Sci. Sta.tist.Cemput.IO 1989 36 52 ) [ 16 ] 張紹良 藤野清次 小柳義夫 全周ノイマン型塊界条件をもっボアソン方程式の数値解法 第 5 回数値流体力学シンポジウム講演論文集 1991 ) 505 508. 堀端康善 正会員 ) 210 川崎市川崎区浮島町 4 1 株 ) 東芝研究開発センター 1979 年東京大学大学院工学系研究科博士課程修了. 工学博士同年 株 ) 東芝入社 現在に至る. 入社後 主として熱と流れの数値シミュレ ション 逆問題の研究に従事 H 本シミュレーション学会 日本数値流体力学会 日本流体力学会 電子情報通信学会 日本気象学会各会員. 1997 年 3 月 25 日受付 ) 一 149 一