将来行ってみたい研究

Transcription

1 中部 CAE 懇話会 : 基礎講座第 7 回 (H 年度 ) 流体伝熱基礎講座 H 年度 4 日目 ( 9 年 8 月 9 日 ( 土 ) ) 流体数値計算法 講師 : 名古屋工業大学 大学院工学研究科教授森西洋平

2 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法

3 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法

4 第 章 時間進行法 目的 : 階常微分方程式に対する時間進行法を理解する. 本章では変数 () に関する次の 階常微分方程式の初期値問題の時間進行法 ( 数値時間積分法 ) を考える. ( ) d d ( ()) ( ) > (.) 階常微分方程式 ( ) (.) 初期条件 つまり 時間刻み幅毎の時間離散点上の値 ( ) ( ) ( ) を 式 (.) を初期条件として 式 (.) を満足させながら求めたい. 4

5 時間進行法を構築するため まず () を () まわりにテイラー展開して式 (.) を使用すると次式を得る. ( ) ( ) d d ( ) d ( ) d () d ( ) d! d L L! d () ( () ) d ( ( ) ) d ( ( ) ) d! d L ( ) ( ( ) ) d L ( ) L! d 時間に関する離散点上の値を表現するため ( ) ( ( ) の表記を使用する. 以降も同様. (.) 5

6 最も単純な時間進行法はテイラー展開の右辺第 項以降を打ち切って構成される陽的オイラー法である. あるいは ( ) d d に対する陽的オイラー法 () () Eler exl 陽的オイラー法では における勾配 と から を計算する. 図. 陽的オイラー法による時間進行 6

7 式 (.) と陽的オイラー法をテイラー展開を通して比べると 打ち切誤差の初項が O( ) となり 陽的オイラー法は 次精度しか持たないことがわかる. 時間進行法の精度を向上させる方法として ) 線形多段階法 ) 段階法 がしばしは採用される. 7

8 ) 線形多段解法 を計算する際に 時刻 および での情報のみならず さらに過去の情報 ( 時刻 - - ) も利用して高次精度化を行う. 計算の出発点近傍で別の時間進行法を使用する必要があり また時間刻み幅の変更が容易でないことに注意が必要. 代表は アダムス バシュフォース (AB) 公式アダムス モルトン (AM) 公式後退差分 (BD) 公式 ) 段解法 : を計算する際に 基本的に時刻 での情報 ( およびそれから計算される値 ) のみを利用して高次精度化を行うもの. 代表は ルンゲ クッタ (RK) 法 以下本章ではこれら時間進行法を取り上げ紹介する. さらに これら時間進行法の数値安定性についても示す. 8

9 補足 : 式 (.) の 階常微分方程式を扱う理由 流体運動の支配方程式は一般に偏微分方程式である輸送方程式として記述されるが なぜ 階常微分方程式を扱うのかを説明する. 例として外力を伴う 次元移流拡散方程式を考える. U x ν x f 変数 (x) が独立変数 と x の双方に依存するので この方程式は偏微分方程式である. またより一般的に 移流速度 UU(x) 拡散係数 νν(x) および外力 ff(x) も時間と空間の双方に依存するものとする. 9

10 ここで 空間に等間隔格子 (x x N) を使用し 離散値を (x ) () のように表現し 空間 階微分を 次精度中心差分で近似する. d d ν f x x U ( N ) 空間離散化後の変数は () なので偏微分が常微分に変わる. さらに変数を [ N ] T とベクトル表記し 端点 N に対する境界条件も含めて上式を書き直すと次の離散システム方程式を得る. d d ( ) () () b() ( () ) A A() は差分の係数を要素とする (N N) の行列 b() は外力と境界条件を含むベクトルである. このように流体運動の支配方程式を空間離散化した後の時間発展方程式は基本的に式 (.) の形となる. 空間多次元の場合も同様である. なお 離散システム方程式の変数はベクトル 式 (.) の変数はスカラーであるが 時間進行法の構成に違いは生じない.

11 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法

12 線形多段階法 d d に対するもの (N) 段の線形多段階法 (Ler ml-se mehods) は一般に次式で記述される. ただし α. α N N β 線形多段階法の一般形 特に α α - β α ( -) とすれば アダムス バシュフォース (Adms-Bshforh) 公式 α α - β α ( -) とすれば アダムス モルトン (Adms-Molo) 公式 α β β ( ) とすれば 後退差分 (Bkwrd-dfferee) 公式

13 アダムス バシュフォース公式 ( ) ( β β β L) 陽解法 (β ) () Adms Bshforh () 図. アダムス バシュフォース公式による時間進行

14 アダムス モルトン公式 ( ) ( β β β β L) 陰解法 (β ) () Adms Molo () 図. アダムス モルトン公式による時間進行 4

15 後退差分公式 ( α α α L) () Bkwrd dfferee 陰解法 (β ) () 図.4 後退差分公式による時間進行 以下でそれぞれに対する 4 次精度までの公式を示す. 5

16 アダムス バシュフォース公式 時間 次精度から 4 次精度までのアダムス バシュフォース (AB) 公式は以下のとおりである. ( ) ( ) ( ) ( ) ( 6 5 ) ( ) ( ) 4 陽的オイラー法 (AB) AB ( アダムス バシュフォース法 ) AB AB4 次精度の公式 (AB) は陽的オイラー法である. 特に精度が明示されていない場合 流体計算では 次精度の公式 (AB) をアダムス バシュフォース法と称する場合が多い. 6

17 7 補足 :AB 公式の導出 (/) ( ) ( ) L β β β とおき 時刻 まわりのテイラー展開 ( ) L 4 6 d d d d d d ( ) ( ) ( ) L L 6 q q d d q d d q d d q から AB 公式が指定した時間精度 ( 次精度ならば O( )) で近似されるように係数 β β - β - を定める. AB 公式を

18 8 補足 :AB 公式の導出 (/) 例として 次精度アダムス バシュフォース法 (AB) を取り上げる. ( ) β β AB の係数 β β - を定めるためテイラー表を利用する. ( ) 6 / β β β β β β d d d d O( ) と O( ) の列それぞれの係数の合計を とする つの条件から 係数 β β - に対する連立方程式が得られる.

19 9 補足 :AB 公式の導出 (/) 係数 β β - に対する連立方程式 β β β これを解き β / β - -/ を得る. これは AB 公式を与える. ( ) ( ) なお テイラー表の 列目 (O( ) ) の合計は AB 公式の打切り誤差の初項 ε を与える. AB ( アダムス バシュフォース法 ) ) ( 5 6 O d d d d β ε

20 補足 : 陽解法による計算 陽解法は既知量の単純代入計算のみで時間進行が行えるので計算コードの作成が容易であり また陰解法に比べて時間ステップあたりの計算量が非常に少なくて済む. 例えば AB は実際には次の形で使用される. ( ) これは既知量 ( 時刻 以前の量 ) の単純代入計算のみで未知量 が得られることを意味し 時間ステップの進行に繰り返し計算を必要としない. このような時間進行法は陽解法と呼ばれる. AB 公式陽的 RK 法

21 アダムス モルトン公式 時間 次精度から 4 次精度までのアダムス モルトン (AM) 公式は以下のとおりである. ( ) ( ) ( ) ( ) ( 5 8 ) 陰的オイラー法 (AM) ( ) ( ) 4 クランク ニコルソン法 (AB) AM AM4 次精度の公式 (AM) は陰的オイラー法であり 次精度の公式 (AM) は通常クランク ニコルソン (Crk-Nolso CN) 法と呼ばれる.AM 公式は の値を得るために ( ) を必要とするが このような時間進行法は陰解法と呼ばれる.

22 補足 :AM 公式の導出 (/) ( ) ( ) L β β β とおき 時刻 まわりのテイラー展開 ( ) L 4 6 d d d d d d ( ) ( ) ( ) L L 6 q q d d q d d q d d q から AM 公式が指定した時間精度 ( 次精度ならば O( )) で近似されるように係数 β β β - を定める. AM 公式を L 6 d d d d d d

23 補足 :AM 公式の導出 (/) 例としてクランク ニコルソン法 (AM) を取り上げる. AM の係数 β β を定めるためテイラー表を利用する. O( ) と O( ) の列それぞれの係数の合計を とする つの条件から 係数 β β に対する連立方程式が得られる. ( ) β β ( ) 6 / β β β β β β d d d d

24 補足 :AM 公式の導出 (/) 係数 β β に対する連立方程式 β β β これを解きβ β /を得る. これはクランク ニコルソン法 (AM) を与える. ( ) ( ) クランク ニコルソン法 (AB) なお テイラー表の 列目 (O( ) ) の合計はクランク ニコルソン法の打切り誤差の初項 ε を与える. d d ε β O( 6 d d ) 4

25 補足 : 陰解法による計算 例えばクランク ニコルソン法 (AM) は実際には上式を変形した次の形で使用される ( ) ( ) 右辺は既知量であるが左辺は未知量 を含む関数であり 上式を満足する を特定するために通常繰り返し計算が必要となる. このような時間進行法は陰解法と呼ばれる. AM 公式後退差分公式陰的 RK 法 5

26 後退差分公式 時間 次精度から 4 次精度までの後退差分 (BD) 公式は以下のとおりである. ( ) ( 4 ) 陰的オイラー法 (BD) BD ( 8 9 ) 6 BD ( ) BD4 次精度の公式 (BD) は陰的オイラー法である. これら後退差分公式はギアー (Ger) 法と呼ばれることもある. BD 公式も陰解法である. 6

27 7 補足 :BD 公式の導出 (/) とおき 時刻 まわりのテイラー展開 から BD 公式が指定した時間精度 ( 次精度ならば O( )) で近似されるように係数 α α α - α - を定める. BD 公式を ( ) L α α α α ( ) ( ) [ ] ( ) [ ] L 6 q d d q d d q d d q ( ) ( ) [ ] ( ) [ ] L 6 q d d q d d q L q

28 8 補足 :BD 公式の導出 (/) 例として 次精度後退差分公式 (BD) を取り上げる. BD の係数 α α α - を定めるためテイラー表を利用する. O( - ) O( ) および O( ) の列それぞれの係数の合計を とする つの条件から 係数 α α α - に対する連立方程式が得られる. ( ) α α α 6 8 / 6 / / d d d d α α α α α α α α α α α α

29 9 補足 :BD 公式の導出 (/) 係数 α α α - に対する連立方程式これを解き α /α -α - / を得る. これは BD 公式を与える. なお テイラー表の 4 列目 (O( ) ) の合計は BD 公式の打切り誤差の初項 ε を与える. BD 4 α α α α α α α ( ) 4 ( ) ) ( 6 8 O d d d d α α ε

30 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法

31 ルンゲ クッタ法 ( ) d d ( ) ( () ) ( > ) (.) (.) 式 (.) を時間進行するための s 段ルンゲ クッタ法 (Rge-K mehod) の一般公式は次式で与えられる. ( ) s s b ただし 係数に対し s ( s) の条件を与える. ( ( ) ) ( s) ( ( ) ) (..) (..) 係数 b および は要求精度を満たすように式 (.) のテイラー展開との比較から決定される.

32 RK 法の係数は次のブッチャー配列 (Bher ble) で表現すると便利である. s ss s s s s s b b b L L M O M M M L L (..) s 段 RK 法のブッチャー配列 T A b

33 なお に対して とするものを ( 陽的 )RK 法と呼ぶ.s 段 ( 陽的 )RK 法の到達可能精度 (s) は次表となることが知られている. s () s () s s 8 9 s 精度と計算量の兼ね合いからは 4 段 4 次精度公式が最も効率が良いことがわかる. また > に対して とするものを半陰的 RK 法と呼ぶ.s 段半陰的 RK 法の到達可能精度は (s)(s) となることが知られている. さらに ある > について のものを陰的 RK 法と呼ぶが s 段陰的 RK 法の到達可能精度は (s)s である.

34 4 補足 : 係数に対する条件式 (..) に対しての必然性はないが このように設定すると扱い易い形となるので 通常この条件を加える. また 常微分方程式の数値解法の教科書における RK 法の一般公式は式 (..) よりも次式によるものが多い. しかしここでは 計算流体力学への応用における明快性から式 (..) を一般公式として採用する. ( ) s s K b s K K

35 5 陽的ルンゲ クッタ法陽的 RK 法は > に対して と設定するものであり s 段陽的 RK 法に対するブッチャー配列は以下となる. 陽的 RK 法のブッチャー配列 s s ss s s s b b b b L L M O O M M M L L L L L L

36 6 この場合 ( として ) ( ) ( ) s ( ) ( ) ( ) ( ) ( ) s と記述すると s 段陽的 RK 法は次のように表現できる ( 更新方式 ). () () () ( ) () ( ) () s b s なおこの公式から は 段目の計算時点での から の間の到達割合いを意味することがわかる.

37 以下 4 次精度までの陽的 RK 法の公式の具体例を挙げる. まず 陽的オイラー法は 段 次精度陽的 RK 法 (RK) と見なすことができ 対応するブッチャー配列と計算式は b () () () となる. これが陽的オイラー法 と一致することは容易に理解できる. 7

38 8 段 次精度陽的 RK 法 (RK) 段陽的 RK 法の一般形は次式で与えられる. () () ( ) ( ) ( ) () ( ) ( ) b b ここで b α を自由パラメータに選ぶと 段 次精度陽的 RK 法 (RK) に対する係数のブッチャー配列は以下で与えられる. α α α α ) /( ) /( b b RK の一般係数 (α は自由パラメータ )

39 特に α の場合を修正オイラー法 α/ の場合を He 法 ( 改良オイラー法 ) と呼ぶ. また α/4 の場合は Rlso の最適公式である. 修正オイラー法 (α) 改良オイラー法 (α/) Rlso の最適公式 (α/4) / / / / / / / 4 / 4 () () RK (α Modfed Eler) () () () () / / 図.5 () () RK( 修正オイラー法 ) による時間進行 9

40 補足 :RK の係数 (/) まず( ()) の時間微分は d d d D D ( D ) d 4 d D 4 d d と書ける. ただし (()) / および D D d ( D ) ( D) ( D) D D とする. これより () のテイラー展開は次式で記述される.! [ ] ( ) ( ) D D ( D ) 4 ( D ) ( D ) ( D ) [ ] D L D 4! 4

41 補足 :RK の係数 (/) 次に RK の一般公式における () () を および () まわりにテイラー展開すると次式を得る. () () ( () ( ) ) L これをRKの時間進行の式 (b () b () ) に代入すると RKは ( ) ( ) ( b b ) ( b b ) L の計算を行っていることに対応する. ただし / である. これを () のテイラー展開と比較して O( ) まで係数を一致させ さらに式 (..) の条件も加えると RK の係数に対する条件式を得る. 4

42 補足 :RK の係数 (/) RK の係数に対する条件式 これは自由度 の連立方程式 ( 式 つ 変数 4 つ ) であるので b α を自由パラメータとして係数の解を表現すると を得る. b b /(α) b -α b / 4

43 4 段 次精度陽的 RK 法 (RK) 段陽的 RK 法の一般形は次式で与えられる. 段陽的 RK 法に対する係数のブッチャー配列は以下となる. 段陽的 RK 法のブッチャー配列 () ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) () ( ) ( ) ( ) b b b b b b

44 代表的な 段 次精度 RK 法 (RK) に対する係数のブッチャー配列を挙げておく. 古典的 RK( / ) HeのRK( /4 /) RlsoのRK( / /4) / 8 /5 / / / 6 8 /5 / 4 / 4 / 5 / / 6 / 4 / 4 / / 4 / 9 / 4 8 / 9 WryのRK( 8/5 /) b となるRK( /6 /4) WllmsoのRK( / /4) / 6 / 4 / 6 7 / 4 / 7 / 4 4 / 7 / / 4 / 4 / / 9 /6 / 6 / 4 / 5 /6 これらのうち b となる RK は更新方式における低容量型 RK で ある. つまり b によって の計算に () が必要なくなるので 更新計算において他の RK 公式よりも配列を つ減らせる.Wry の RK は積み上げ方式における低容量型 RK である.Wllmso の RK は Wllmso の計算法に従えば低容量 (N 容量 ) 型となる RK である. / / / 4 / 9 8 /5 44

45 補足 :RK の係数 (/) RK と同様にして RK に対する係数の条件式も得られる. b b b b b b b / b / 6 / これは自由度 の連立方程式 ( 式 6 つ 変数 8 つ ) である. 45

46 46 補足 :RK の係数 (/) を得る. また b α を自由パラメータとして 以下の解が知られている. と を自由パラメータとして係数の解を表現すると / に対して ( ) ( ) 6 b ( ) 6 b ( ) 6 b α α α α 4 / 4 / ) /(4 ) /(4 / / / / α α α α 4 / 4 / ) /(4 ) /(4 / /

47 47 4 段 4 次精度陽的 RK 法 (RK4) 4 段陽的 RK 法の一般形は次式で与えられる. 4 段陽的 RK 法に対する係数のブッチャー配列は以下となる. 4 段陽的 RK 法のブッチャー配列 () ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) b b b b b b b b

48 4 次精度を満足する 4 段陽的 RK 法の係数の一般解はいくつか知られているが ここでは代表的な 4 段 4 次精度 RK 法 (RK4) に対する係数のブッチャー配列を挙げておく. 古典的 RK4 K の RK4(K の /8 公式 ) / / / / 6 / / / / 6 / / / / /8 /8 /8 /8 Gll の RK4 / / / ( )/ ( ) / / ( ) / / 6 ( ) / 6 ( ) / 6 / 6 48

49 補足 :RK4の係数(/) RKあるいはRKと同様にしてRK4に対する係数の条件式も得られる b b b b b b b b b b b4 b b44 b b44 b b44 / / 4 / ( ) / 6 b4 4 4 ( ) / b4 4 4 ( ) / 8 b / 4 これは自由度 の連立方程式 ( 式 変数 ) である. 49

50 5 補足 :RK4 の係数 (/) を得る. この解は / / で K の RK4 を与える. まず と を自由パラメータとし ( -)( -)( - ) の場合について RK4 の係数の解を表現すると 4 ( ) ( ) 4 ( ) ( ) ( ) ( ) [ ] ( ) ( ) [ ] ( )( )( ) ( ) ( ) [ ] ( ) b ( )( ) b ( )( ) b ( ) ( ) 4 b b b b また

51 補足 :RK4 の係数 (/) /に対してはb α を自由パラメータとする次のパラメータ 解が知られている. / / / / /(6α ) /(6α ) α α / 6 / α α / 6 この解は α/ で古典的 RK4α( )/6 で Gll の RK4 を与える. 他にも 4 つの パラメータ解が知られているが そのうちの つを示しておく ( ただし α ). / / / / 8 / 8 /(α ) /(α ) /(4α ) /(α ) /(α ) / 6α / 6α / 6 / 6 α / α / 6 α / α / 6 5

52 陰的ルンゲ クッタ法 陰的 RK 法はある>について なので 基本的に式 (..) をそのまま適用する必要がある. 例えば 段陰的 R K 法に対するブッチャー配列と計算式は および b () () b { ( () ) ( () )} { ( () ) ( () )} { ( () ) ( () )} b b であり () と () の連立方程式を解いた後に を計算する. 一般にs 段陰的 RK 法では () () (s) の連立方程式をニュートン法等の繰り返し法で時間ステップ毎に解く必要があり計算量が膨大となる. 5

53 半陰的 RK 法は>について と設定するものであり 例えば 段半陰的 RK 法に対するブッチャー配列と計算式は および () () b b () ( ) () () { ( ) ( )} { ( () ) ( () )} b b となる. この場合は () および () に対する単独の陰的方程式を順番に解いた後に を計算すれば良いので 段あたりの計算負荷は AM 法と同程度となる. 陰的 RK 法はs 段で最大 s 次精度まで 半陰的 RK 法ではs 段で最大 (s) 次精度まで到達可能である. 5

54 陰的 RK 法 (IRK 法 ) および半陰的 RK 法 (SIRK 法 ) の例を挙げておく. 陰的オイラー法 ( 段 次 IRK 法 ) 陰的中点則 ( 段 次精度 IRK 法 ) CN 法 ( 段 次 SIRK 法 ) / / / / / / 段 4 次 IRK 法 段 次 SIRK 法 ( ノルセット法 ) ( ) / 6 / 4 ( ) / ( ) / 6 ( ) / 6 ( ) / 6 ( ) / / 4 ( ) / 6 / ( ) / / / / / 6 段 6 次 IRK 法 ( 5 5) / ( 5 5) / / 5/6 5/6 5 / 4 5/6 5 / 5/8 /9 5 /5 /9 /9 5 /5 4/9 5/6 5 / 5/6 5 / 4 5/6 5/8 54

55 ここで 陰的オイラー法 陰的中点法 クランク ニコルソン法 (CN 法 ) について補足しておく. 陰的オイラー法は 段 次精度陰的 RK 法とみなせ ブッチャー配列と計算式はそれぞれ および と記述できる. この場合の つの計算式は等価で () であり これを用いると AM 公式での陰的オイラー法 () が得られる. () ( ) () ( ) 55

56 56 陰的中点則は 段 次精度陰的 RK 法とみなせ ブッチャー配列と計算式はそれぞれ と記述できる. つの計算式を比べると () ( )/ となっており これが陰的中点則 であることがわかる. および / / () () ) (

57 クランク ニコルソン法 (CN 法 ) は 段 次精度陰的 RK 法とみなせ ブッチャー配列と計算式はそれぞれ および / / () () () () ( ) ( ) () () ( ) ( ) と記述できる. この場合 () () であり これが AM 公式でのクランク ニコルソン法 であることがわかる. ( ) 57

58 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法 58

59 時間進行法の安定性解析 線形テスト問題 d d λ ( ) λ λ R λ I に対する時間進行法の数値安定性を考える. まずこの問題の解析解は λ λ ( ) e e R λ e I で与えられ λ R ならば > で () は有界となる. λ I Sble Usble λ R 図.6 線形テスト問題の安定領域 ( 複素平面 ) 59

60 補足 : なぜ線形テスト問題を考えるのか (/) 流体運動の支配方程式は輸送方程式であるので その本質を表現するモデル方程式 ( 次元の移流 拡散方程式の空間離散式 ) を考える. d d δ δ U ν δx δx ここで U: 移流速度 ν: 拡散係数 δ/δx: / x に対する空間離散オペレータ 解析的な扱いを容易にするため x 方向に周期的な場を仮定し 離散フーリエ変換を導入する. ˆ k N π () ( x ) ex[ ω ( k) x ] k N / ~N / ω( k) k L N x N / ( x ) ˆ () ex[ ω( k ) x ] ~N k k N / L x x 6

61 補足 : なぜ線形テスト問題を考えるのか (/) 離散フーリエ変換によりモデル方程式を対角化する. dˆ k ( ) [ νω ( k) Uω ( k) ] ˆ k ( k) k N / ~N / d ここで ω (k) とω (k) は空間 階および 階微分に対する離散オペレータの修正波数で 空間離散化手法に応じて定まる. 例えば ( k ) ω( k ) ω ( k ) ω( k ) フーリエ スペクトル法 : ω 次精度中心差分 : ω ( k ) s ( ω( k ) x) x ω ( k ) s ( ω( k ) x / ) ( x / ) 以上より モデル方程式は線形テスト問題 ( 固有値問題 ) に帰着され 固有値の実部は拡散項に 固有値の虚部は移流項に関係することが理解される. dˆ k d ( ) R I R I ( ) λ λ λ λ ν ω ( k ) λ U ω ( k ) λ ˆ k k k k k なお 純粋移流問題の安定性はω (k) の最大値 純粋拡散問題の安定性はω (k) の最大値に対応して定まる. k k 6

62 線形テスト問題の解の有界性が数値計算においてどのような条件で満足されるかを検討するため まず陽的オイラー法を線形テスト問題に適用してみる. ( ) λ ( λ) この操作は 初期値 数値解の複素増幅率を ξ として ξ と記述するとき 陽的オイラー法の複素増幅率が ( z) ξ ただし z λ となることを意味する. ここで zλλ R λ I としている. 6

63 これから明らかなとおり 陽的オイラー法による数値解が有界であるための条件は ξ z で与えられる. これを満足するzは複素平面において (-) を中心とする半径 の円の内側に存在する. この領域を陽的オイラー法の絶対安定領域と呼ぶ. これは 陽的オイラー法を用いて安定な数値解を得るためには 時間刻み幅の設定に制限があることを示している. 具体的には λを負の実数とすると λ の範囲に時間刻み幅 を設定する必要がある. λ I Eler exl Usble Sble λ R 図.7 陽的オイラー法の絶対安定領域 6

64 次に 陰的オイラー法を線形テスト問題に適用し絶対安定領域を調べる. ( ) λ これより 陰的オイラー法の複素増幅率 ξ は次式で与えられる. ξ ( z) ただし z λ これに対応する絶対安定領域は次式で与えられる. ξ z これを満足する z は複素平面において (-) を中心とする半径 の円の外側の全領域である. この領域は λ R となる複素平面の左半分を全て含むので 陰的オイラー法は λ R に対し時間刻み幅 によらず無条件安定な数値解を与える. λ I Sble λ R Eler ml Usble 図.8 陰的オイラー法の絶対安定領域 64

65 陰的オイラー法のように絶対安定領域が複素平面の左半分 (λ R ) 全てを含む時間進行法は A 安定 ( 絶対安定 ) という. Sble Usble なお 絶対安定の緩和表現として 絶対安定領域が複素平面上で無限扇形領域 {z; -α < (π-rg z) < α} (<α<π/) を含む場合 A(α) 安定と呼ぶ. α α Im(z) O Re(z) 図.9 A(α) 安定 ところで 陽的オイラー法および陰的オイラー法は 次精度の時間進行法であるため打ち切り誤差の影響が懸念される 65

66 66 そこで 次精度の時間進行法である中点蛙跳び法 (le-frog 法 ) を線形テスト問題に適用してみる. これに ξ および zλ を代入して整理すると 中点蛙跳び法の複素増幅率 ξ はの解として次式で得られる. ここでは複素増幅率 ξ の解を求めるのではなく 絶対安定領域の境界 ξ を与える z の範囲を調べてみる. そのため上式を変形し ξe θ を代入すると ( ) λ ξ ξ z ( ) θ ξ ξ θ θ θ θ os e e e e z を得る. これは z が虚軸上の (- ) の範囲のみで安定 つまり中点蛙跳び法は λ R 以外では数値解が無条件に不安定となることを示している. 以上のように 時間進行法は必ずしも精度が高いものが適切なわけではなく 数値安定性も含めて検討する必要がある.

67 線形多段階法の絶対安定領域 線形多段階法の一般形 ( ただし α ) α N N β を線形テスト問題に適用し ξ および zλ を代入すると次式を得る. α ξ z β ξ N N この式による複素増幅率の根が全て ξ を満足する場合に線形多段階法は絶対安定である. 67

68 線形多段階法の安定限界曲線 ( 絶対安定領域の境界 ) つまり ξ を与える z は LM 法の安定限界関数 z ϕ ( ξ ) N N α β ξ ξ LM 法の安定限界関数 に対して ξe θ を代入し から π の範囲の θ に対する z を描けば良い. 以下に AM 公式 AM 公式および BD 公式に対する安定限界関数と安定限界曲線を示しておく. 68

69 AB 公式 z ξ 陽的オイラー法 (AB) z z z ( ξ ) ξ ξ ξ ξ 55ξ ( ξ ) 6ξ 5 4ξ 59ξ ( ξ ) AB 7ξ 9 AB AB4 λ I AB Eler exl s AB s AB4 s s λ R 図. AB 法の安定限界曲線 69

70 AM 公式 z z z z ξ ξ ( ) ξ ξ ξ 5ξ 9ξ ( ξ ) 8ξ 4ξ 9ξ 陰的オイラー法 (AM) クランク ニコルソン法 (AM) ( ξ ) AM 5ξ AM4 λ I 4 s AM AM4 s AM s s AM s Eler exl λ R 図. AM 法の安定限界曲線 7

71 BD 公式 ξ z ξ z z z ξ 4ξ ξ 4 ξ 8ξ 6ξ 陰的オイラー法 (BD) BD 9ξ ξ 4 BD 5ξ 48ξ 6ξ 6ξ 6 λ I 7 o 86 o BD Eler exl s BD s s s λ R s BD BD4 図. BD 法の安定限界曲線 BD4 λ I BD BD4 s s λ R 図. BD と BD4 の原点付近の安定限界曲線 7

72 以上より 線形多段解法 (AB 法 AM 法 BD 法 ) では 精度が上がると絶対安定領域が狭くなることがわかる. また これらの中で A 安定な時間進行法は 陰的オイラー法 (AM) クランク ニコルソン法 (AM) および 次精度後退差分法 (BD) のみである. なお BD は A(86 ) 安定 BD4 は A(7 ) 安定である. 7

73 7 RK 法の絶対安定領域線形テスト問題に s 段 RK 法を適用し行列表示すると次式を得る. ただし zλ である. s s s ss ss s s s s s s s s s zb zb z z z z z z z z z z z z z z M M L L L L L L M M O O O M L L L ) ( ) ( () () ) ( ) ( ) ( ) (

74 この連立方程式を Crmer の公式で解き ( 解の s 成分 ) を求めると次式を得る. ( T ) I za z de R( z) R( z) eb de ( I za) ここで s s 行列 A と I の成分を (A) および (I) δ また s 元ベクトル b と e の成分を b T (b b s ) および e T ( ) としている. R(z) は RK 法の安定性因子と呼ばれ RK 法の絶対安定領域は R(z) で与えられる. 74

75 補足 :RK 法の安定性因子の導出 (/) まず RK 法に対するブッチャー配列を A b T と表現し s s 行列 A と I の成分を (A) トルの成分を b T (b b s ) T ( s ) e T ( ) T ( ) T ( () (s) ) および (I) δ また s 元ベク とすると 線形テスト問題に適用された s 段 RK 法は次のように記述される. ( I za) T zb e 75

76 補足 :RK 法の安定性因子の導出 (/) 上記連立方程式をCrmerの公式で解いた の解 ( 解のs 成分 ) はまず次式となる. ( I za) de T zb ( I za) de T zb e ここで ブロック行列の行列式に関する数学公式 A C ( ) ( ) de de D de A CD B B D を適用すると 分母および分子に現れている行列式はそれぞれ ( I za) de T zb de( I za) となるので ( T de I za z ) de( I za) を得る. ( I za) e de T zb de R( z) T ( I za zeb ) ( T ) eb de I za zeb R( z) de ( I za) 76

77 陽的 RK 法に対しては >I で なので de(i-za) であり 安定性因子 R(z) は ( T I za zeb ) R( z) de 陽的 RK 法の安定性因子 与えられる. これは一般に z の s 次多項式 R ( z) q z q z L q s z となり 陽的 RK 法ではA 安定なものが存在しないこと ( z - で R(z) となるため ) を示している. また 線形テスト問題に対する の解析解のテイラー展開 λ e z z! s m z L m! L と比べると s 段 次精度陽的 RK 法の安定性因子は z z R ( z) z L q z L!! となることが解る. q s z s 77

78 特に s 4 に対しては s 段 s 次精度陽的 RK 法が構成できるが この場合の安定性因子はそれぞれ以下となる R ( z) z R ( z) z z! z R ( z) z! z! z z R ( z) z!! 陽的オイラー法 (RK) RK z RK 4 4! RK4 陽的 RK 法に対する安定限界曲線を図.4 に示す. 線形多段階法と異なり 陽的 RK 法は精度が上がると絶対安定領域が広がる. λ I RK RK4 s s RK 4 s s λ R Eler exl 図.4 陽的 RK 法の安定限界曲線 78

79 79 陰的および半陰的 RK 法の安定性因子は z の有理関数となるので A 安定となる可能性がある. 節で紹介した陰的および半陰的 RK 法は全て A 安定であり それぞれに対する安定性因子は以下となる. z z R ) ( ) ( z z z R 6 6 ) ( z z z z R ) ( z z z z z R ) ( z z z z z z z R 陰的オイラー法 ( 段 次 IRK) 陰的中点法 ( 段 次 IRK 法 ) および CN 法 ( 段 次 SIRK 法 ) 段 次 SIRK 法 ( ノルセット法 ) 段 4 次 IRK 法 段 6 次 IRK 法

80 輸送法方程式の時間進行を考える場合 移流項の数値安定性は固有値の虚部に また拡散項の数値安定性は固有値の実部に支配されるので 時間進行法の安定性のまとめとして 実数および純虚数の固有値に対する安定領域を表. に示しておく. 8

81 表. 実数および純虚数の固有値に対する安定領域 (ABAMBDRKIRK) AB(Eler exl) AB AB AB4 AM(Eler ml) AM(Crk - Nolso) AM AM4 BD BD BD4 RK RK RK4 Rel ( λ R - λ R - λ λ -. λ Sble Sble Sble Sble R R R -6 λ R - λ Sble Sble Sble R - λ R -.5 λ R -.79 λ ) Imgry R Usble( λ ) I Usble( λ ) λ λ.7.4 Usble( λ ) I Usble( λ ) Sble I λ.94 I λ 4.7 I Usble( λ ) λ λ.7.8 Remrks A - sble A - sble A - sble A(86 A(7 - sge IRK Sble Sble A - sble - sge IRK4 Sble Sble A - sble - sge IRK6 Sble Sble A - sble I I I I I I o o ) - sble ) - sble 8

82 なお 時間進行法の安定性に関して次の定理が知られている. ダルキスト (Dhlqs) の定理 : ) 陽的ルンゲ クッタ法は A 安定にはならない. ) 陽的な線形多段解法は A 安定にはならない. ) A 安定な線形多段階法の精度は を超えない. 4) 次精度で A 安定な線形多段解法のうち 局所誤差の絶対値が最も小さいものはクランク ニコルソン法 (AM) である. 8

83 第 章の参考文献 () 三井斌友 数値解析入門 常微分方程式を中心に - 朝倉書店 (985). () 三井斌友 常微分方程式の数値解法 岩波書店 (). () U.M. アッシャー L.R. ペツォオルド共著 中森眞理尾雄監修 嘉村友作 三ツ間均共訳 常微分方程式と微分代数方程式の数値解法 培風館 (6). (4) Ger C.W. Nmerl Il Vle Problems Ordry Dfferel Eqos Pree-Hll (97). (5) l S.O. Nmerl Mehods for Il Vle Problems Ordry Dfferel Eqos Adem Press (988). (6) Bher J.C. Nmerl Mehods for Ordry Dfferel Eqos Joh Wley d Sos (). (7) Wllmso J.H. Low-Sorge Rge-K Shemes J. Com. Phys. Vol.5 (98)

84 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法 84

85 第 章 非圧縮性流れの計算アルゴリズム 目的 : 非圧縮性流れの計算アルゴリズム ( 圧力ベース解法 ) を理解する. 非圧縮性流れの運動はナビエ ストークス (NS) 方程式および連続の式を用いて記述される. v v ( vv) ν v s (.) NS 方程式 (.) 連続の式 上式中で は時間 v は速度ベクトル は圧力 / 密度 ( 本章では以降圧力と呼ぶ )ν は動粘度 は勾配演算子 s は外力等によるソース項のベクトルである. 動粘度は定数とする. また 式 (.) の左辺第 項は移流項 右辺第 項は圧力勾配項 右辺第 項は粘性項と呼ばれる. 85

86 本章では 式 (.) および式 (.) を解いて速度 v および圧力 の解を求める非圧縮性流れの一般的な解法を示す. 原始関数法 空間および時間離散化手法を定めれば式 (.) および式 (.) は離散化できるが 非圧縮性流れでは連続の式 (.) が NS 方程式 (.) に対する拘束条件となるので これらの式が同時に満足されるように速度 v および圧力 の数値解を求める解法が必要となる. 86

87 速度 v については式 (.) を基に時間発展を行えば良いが 非圧縮性流れでは圧力 に関する時間発展方程式が存在しない. このため 通常は式 (.) に発散 ( ) を作用させて得られる圧力のポアソン方程式 ( v) [ v ( vv) s] ν (.) 圧力の ポアソン方程式 あるいはこれに関連した方程式と連続の式 (.) とを基に圧力を求めながら式 (.) の時間発展を行い速度 v が求められる. 非圧縮性流れの計算アルゴリズム ( 圧力ベース解法 ) 87

88 NS 方程式 (.) の時間進行法に応じてそれぞれ経済的な解法がこれまでに開発されている. 表. 非圧縮性流れの計算アルゴリズム 計算アルゴリズム NS 式の各項の時間進行法 移流項粘性項圧力項 MAC 法系統陽解法陽解法 フラクショナル ステップ法系統陽解法陰解法 陰解法陽解法 ( 陰解法 ) SIMPLE 法系統陰解法陰解法 それぞれの計算アルゴリズムについて 節 節 節で説明する. 88

89 89 解法の説明に具体的な空間離散式が必要となる場合には 次元の流れ場を考える. その場合の NS 方程式 (.) および連続の式 (.) は次式で与えられる. 上式中で v は速度ベクトル v の xy 方向成分 s x x y はソース項ベクトル s の xy 方向成分である. s x y x x y v x ν s y y v x v y y vv x v v ν y v x (.4) (.4b) (.5) 次元流れの基礎式

90 離散値の表現方法として (x y ) 点上に配置された変数の離散値を例えば と表記する. さらに 時間離散点 における を と表記する. は時間刻み幅である. 節で示されるとおり非圧縮性流体の差分解析ではスタガード格子の使用が推奨されるので 本章での各解法の説明は基本的にスタガード格子を用いて行われる. レギュラー格子スタガード格子コロケート格子 9

91 なお 有限体積法に基づいて空間離散化が行なわれる場合には 通常式 (.) の代りに次式 v { [ ] T } s ( vv) ν ( v) ( v) (.6) が基礎式として用いられる.( v) T は速度勾配テンソル ( v) の転置テンソルである. 粘性係数 νが定数の場合には連続の式 (.) から式 (.) と式 (.6) とは等しくなる. 渦動粘度 ν を用いる乱流モデルの運動方程式は式 (.6) のνを (νν ) で置き換えたものになり またνが空間的な変数である場合も扱えるので 非圧縮性流れに対する運動方程式としては式 (.) よりも式 (.6) の方が一般的である. しかし式 (.6) を用いると解法の統一的な説明が煩雑になるので 節 節 節の解法の説明では式 (.) を運動方程式として用いる. 乱流モデルを用いる場合および有限体積法により空間的離散化が行われる場合に対しては 粘性項の空間離散式が式 (.6) 中の形を基に構成されているものと解釈して頂きたい. 9

92 また 輸送量 φ の輸送方程式の一般形は次式となる. φ ( vφ ) ( Γφ φ ) sφ (.7) ここでΓ φ はφの拡散係数 s φ はソース項である. 式 (.7) の形の輸送方程式を流れ場と同時に解く場合 通常時間進行の各ステップにおいて速度 vおよび圧力 が求められた後に式 (.7) を単純に時間進行して φが求められるので 式 (.7) の時間進行法については特に説明しない. 9

93 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法 9

94 差分格子系と MAC 法系統の解法 NS 方程式 (.) の移流項と粘性項を陽的に時間離散化し 連続の式 (.) を新たな時刻 において満足させると次式を得る. v v ( v v ) ν v s (.8) v (.9) 式 (.8) では時間に関して 次精度の陽解法が用いられているが MAC 法系統の解法で時間精度の向上を図りたい場合には移流項および粘性項に対してより高次精度の陽的時間進行法 ( 例えば 次精度アダムス バシュフォース法 ) が用いられる. 94

95 式 (.8) に発散 ( ) を作用させて式 (.9) を用いると圧力 のポアソン方程式が得られる. [ v ( v v ) s] v ν (.) MAC 法系統の解法では 基本的には圧力のポアソン方程式 (.) を解いて式 (.9) を満足する速度場 v を式 (.8) から求める. 式 (.) を解いて を得る 式 (.8) から v を得る この流れ場は式 (.9) を満足する 95

96 レギュラー格子を用いる差分法の問題点 速度成分 v および圧力 を同じ空間離散点上に配置する差分法は基礎方程式の離散化が最も容易である. この空間配置による差分格子をレギュラー格子と呼ぶ ( 図.). y y : v y y y y x x x x x x 図. レギュラー格子 96

97 97 式 (.8) および式 (.9) を 次元で記述し 空間微分をレギュラー格子上の差分で近似して整理すると次式を得る. f x x g y y v v y y v v x x (.) (.b) (.) ここで f および g は s x y x y v x f ν s y y v x v y vv x v g ν であり f と g はそれぞれ (x y ) 点上での f と g の離散近似式である. x y x x x x x y y y y y : v

98 98 圧力 は 式 (.) により および v が時間進行された後に式 (.) を満足するものでなくてはならない. そのような圧力は次式を満足する. ( ) R y y g g x x f f D y y y y y y x x x x x x (.) x y x x x x x y y y y y : : v : 図. レギュラー格子の問題点ここで ( ). R y y v v x x D 式 (.) は式 (.) を式 (.) に代入して得られる. しかしながら式 (.) の左辺は隣り合う格子点の値を参照しない離散式であるため ( 図.) この式から得られる圧力の数値解は空間的に つ飛びの振動を起こす可能性があり滑らかな物理的な圧力解が得られない場合がある.

99 99 式 (.) は式 (.) を空間的に離散化したものと見なせるので 式 (.) の左辺を 階微分に対する次の差分式で置き換えて圧力を求めることも考えられる. (.4) 上式による圧力 の解は滑らかで物理的なものとなるが 逆にこの圧力を用いて計算される速度 および v は連続の式の離散式 (.) を満足しなくなる. つまり レギュラ格子を用いた差分法による非圧縮性流体の解法では 圧力の解の物理性あるいは連続の式の満足度のいずれかが犠牲になる. ( ) R y y g g x x f f D y y y y y y x x x x x x

100 MAC 法とスタガード格子 レギュラー格子の問題点を解消する差分格子として 非圧縮性流体の数値解析ではスタガード格子が用いられる ( 図.). x / スタガード格子では 圧力の配置点に対しそれぞれの方向に半格子ずれた位置に速度成分が配置される. つまり 圧力 の配置点を (x y ) とすると x 方向速度成分 は点 (x / y )y 方向速度成分 v は点 (x y / ) 上に配置される. y y y y y x x / x x / x x x / x 図. スタガード格子 x / y / y / y / y / : : v : y /

101 x / x / x / x / x / y / y y / y / NS 方程式 (.8) の各方向成分はそれぞれの速度成分の配置点上で また連続の式 (.9) は圧力の配置点で空間的に離散化される. y y y y x x x x x y / y / : : v : / / x / f / (.5) v / v / g / y / (.5b) / x / v / v y / (.6)

102 式 (.5) を式 (.6) に代入すると次式が得られる. ( ) S y g g x f f D y y y x x x / / / / / / / / ここで ( ) S y v v x D / / / / (.7) 式 (.7) の左辺は式 (.4) の左辺と等価であり これを解いて求められる圧力の数値解 は滑らかで物理的であり さらにこの圧力を用いて式 (.5) により時間発展される速度 および v は式 (.6) を離散的に満足する.

103 スタガード格子の使用が提案されたこの解法は MAC 法 (Mrker Ad Cell 法 ) と呼ばれる (). 原論文においてはマーカーを飛ばして自由表面の位置を決定しながら流れ場を解く解法として MAC 法が提案されているが 現在では流れ場の解法の部分のみに対して MAC 法の名称が用いられている. MAC 法の計算手順 ) 初期条件として と v を与え ) 式 (.7) を解いて を求めた後 ) 式 (.5) から と v を求め 4) と v を と v として ) に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す.

104 SMAC 法 MAC 法では圧力のポアソン方程式 (..7) の右辺の境界値および圧力の境界条件を適切に与えることが一般に煩雑である. これを解消する解法にSMAC 法 (Smlfed MAC 法 ) がある ().SMAC 法による解法を説明するため まず式 (.5) を次のように分解する. 中間速度 ˆ / / f / x / ˆ / v / g / y / v (.8) (.8b) / ˆ / δ x δ / (.9) また v / v ˆ / δ δ y δ / (.) (.9b) δ は圧力補正量 4

105 5 式 (.9) を式 (.6) に代入すると δ のポアソン方程式が得られる. ( ) S D y y y x x x / / / / ˆ δ δ δ δ δ δ δ δ ここで ( ) S y v v x D / / / / ˆ ˆ ˆ ˆ ˆ (.)

106 SMAC 法の計算手順 ) 初期条件として v を与え ) 式 (.8) から^ と^v を計算し ) 式 (.) を解いてδを求めた後 4) 式 (.9) から とv および式 (.) から を求め 5) v を v として ) に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. SMAC 法 ( および MAC 法 ) ではポアソン方程式を解くために計算時間の大部分が費やされるので 楕円型方程式の離散式が経済的に解ける行列解法を使用する必要がある. 6

107 MAC 法と比べた場合の SMAC 法の利点は境界条件の設定がわかり易い点にある. 一例として 境界 (x N/ y ) での速度の値 BC N/ が与えられている場合を考える ( 図.4). この場合 まず中間速度 ^ に対する境界値として BC ˆ N / N / を与える. (.) 境界条件は N/ BC N/ であることを要求するので 式 (.9) から δ N δ N (.) が与えられる. y x x N / x N x N/ x N/ BC N/ y δ N δ N 図.4 SMAC 法での境界条件 7

108 補足 : 行列分解による SMAC 法の計算アルゴリズムの導出 (/) 速度 (v) 圧力 および (fg) の離散値ベクトルをそれぞれ および f と表記すると 式 (.5) および式 (.6) は次式のように表現される. D b Gδ G f ただし 圧力は δ としている. また G は離散的勾配演算子 D は離散的発散演算子であり これら離散的空間演算子では境界条件も考慮されている. 境界条件から生じる項は NS 方程式および連続の式それぞれに対して b および b に含まれているものとする. b これらをまとめて行列表示すれば次式となる. G D δ I ( G f b ) b 8

109 補足 : 行列分解による SMAC 法の計算アルゴリズムの導出 (/) ここで 行列表示された基礎式の左辺の行列は次のように LU 分解できる. D I DG G I δ I これはさらに次の 式へと分解できる. ( G f b ) b ( G f b ) I ˆ D DG δˆ b I G ˆ I δ δˆ これらの行列を展開して整理すると SMAC 法の計算アルゴリズムを得る. ( G f ) ˆ b DGδ ( D ˆ b ) DG : 離散ラプラス演算子 ˆ Gδ SMAC 法では分離誤差は生じない 9

110 4 HSMAC 法 ( または SOLA 法 ) MAC 法および SMAC 法の計算時間の大部分はポアソン方程式の離散式を解くために費やされる.HSMAC 法 (Hghly Smlfed MAC 法 )( または SOLA 法 ) () は δ のポアソン方程式 (.) を解かずに速度と圧力の同時補正計算により連続の式を満足させる SMAC 法の近似解法である. まず 式 (.) の左辺の対角項 ( δ に関する項 ) のみを残して δ の式を作り ( 優対角近似 ) 緩和係数 β を乗じると次式を得る. δ x x x β ( Dˆ ) S y y y / / / / (.4)

111 式 (.9) および式 (.) から理解されるとおり 点 () 上の圧力補正量 δ のみが与えられた場合 点 () まわりの速度 / -/ v / v -/ および圧力 のみがその影響を受ける. これより 点 () まわりの速度および圧力をそれぞれ次式により補正する. ˆ / ˆ / δ x / ˆ / ˆ / δ x/ v ˆ / vˆ / δ y / v ˆ / vˆ / δ y / および δ (.5) (.5b) (.5) (.5d) (.6) y y y y y x x / x x / x x x / x / x x / y / y / y / y / : : v : y /

112 HSMAC 法の計算手順 ) 初期条件として v を与え ) 式 (.8) から^ と^v を計算し ) 式 (.4) 式 (.5) 式 (.6) による点 () まわりの速度および圧力の補正計算を全格子点に対して実行し さらにそれをMAX (D^ S ) <εが満足されるまで繰り返す (εは連続の式に対する許容誤差). MAX (D^ S ) <εが満足された時点での速度と圧力を v とする. 4) v を v として ) に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. なお式 (.4) から 連続の式 (D^ S ) が満足されれば速度と圧力は補正されないことがわかる.

113 MAC 法および SMAC 法と比べた場合の HSMAC 法の利点は圧力のポアソン方程式を解く必要が無いという点にある. しかし HSMAC 法は SMAC 法の近似解法であり解の収束が式 (.4) 中の緩和係数 β の値に大きく依存することが欠点となる. なお β には <β の範囲の値が用いられ 経験的に β.7 程度の値が推奨されている ().

114 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法 4

115 フラクショナル ステップ法 先の MAC 法系統の解法では NS 方程式の粘性項および移流項が陽的に時間進行されているので数値安定性条件から時間刻み幅に対する制限が生じる. 特に壁面近傍に差分格子を集中させて速度の滑り無し条件を課す場合 MAC 法系統の解法では粘性項に対する数値安定性条件が厳しくなり非常に小さな時間刻み幅を使用する必要が生じる. 数値安定性条件の問題を解決するには陰的な時間進行法を導入すれば良いので ここでは NS 方程式の粘性項に対してクランク ニコルソン法を適用する解法を考える. 5

116 この場合 式 (.) と式 (.) を時間離散化した方程式は次式となる. v v ν ( v v ) ( v v ) v v s v ν (.7) (.8) NS 方程式の移流項には 次精度のアダムス バシュフォース法が適用されている. ところで 式 (.7) および式 (.8) を直接解いて速度 v および圧力 の連立解を得るには大規模な行列計算が必要となるので経済的ではない. これは NS 方程式の粘性項あるいは移流項を陰解に時間進行する解法に共通の問題点である. フラクショナル ステップ法は 粘性項と圧力項の分離に基づき非圧縮性流れの陰的な非定常数値計算を経済的に実施可能とする解法である. 6

117 フラクショナル ステップ法と分離誤差 まず フラクショナル ステップ法の一例として Km & Mo (KM) () によるフラクショナル ステップ法を示す. この解法では 式 (.7) は次の様に分離される. vˆ v v vˆ ϕ ( ) ( v v v v ) vˆ v s ν ν (.9) (.) 式 (.7) と式 (.9) との差をとり式 (.) と比べるとポテンシャル関数 ϕ と圧力 との関係が ν ϕ ϕ (.) で与えられることがわかる. 式 (.) と式 (.8) からポテンシャル関数 ϕ のポアソン方程式が得られる. ϕ vˆ (.) 7

118 KMのフラクショナル ステップ法の計算手順 ) 初期条件としてv を与え ) 式 (.9) から^v を計算し ) 式 (.) を解いてϕ を求めた後 4) 式 (.) からv を求め 圧力が必要ならば式 (.) から 求め 5)v をv として ) に戻る. 以上の手順を必要な時間まで繰り返す. ところで KMのフラクショナル ステップ法では式 (.7) を式 (.9) および式 (.) へと分離して解法が構成されているが このような分離の際には分離誤差が生じてしまう. 加えて 上記のように空間的に離散化されていない方程式を基にフラクショナル ステップ法を構成すると中間速度 ^v の境界条件の設定が曖昧となってしまう. 8

119 分離誤差および中間速度の境界条件設定の問題を解決するため 時間および空間に関して完全に離散化された方程式を基に解法を議論する. ここでは式 (.7) と式 (.8) の空間微分演算子についても離散化した方程式を考える. ν G L ν L D b s b (.) (.4) ただし 速度 (v) および圧力 の離散値ベクトルをそれぞれ および 移流項の離散ベクトルをと表記している. また Lは離散的ラプラス演算子 Gは離散的勾配演算子Dは離散的発散演算子である. これらの離散演算子は境界条件も考慮されたものであり 境界条件から生じる陽的な項はN S 方程式および連続の式それぞれに対してb およびb で表現されている. 9

120 式 (.) と式 (.4) を行列を用いて表示すれば次式となる. KM KM D G A r r (.5) ここで A は A(I-.5νL) I は単位行列であり また式 (.5) の右辺のベクトルは b b s r r KM KM L I ν で定義される. 式 (.5) の右辺は陽的に扱う項および境界条件により生じる項 左辺が陰的に扱う項である.

121 フラクショナル ステップ法の導入に伴う分離誤差の評価を行うため 式 (.5) に対して次の近似式を導入する (). (.6) 上式中の行列 B は A - に対する近似行列である. 式 (.6) の左辺の行列は次式のように LU 分解できる. ( ) KM KM D G AB A r r KM KM I BG I DBG D A r r (.7) 上式はさらに次の 式へと分離できる. KM KM DBG D A r r ˆ ˆ ˆ ˆ I BG I (.8) (.9)

122 式 (.8) および式 (.9) の行列を展開して整理すると次式を得る. A ˆ r KM (.4) DBG ( Dˆ b ) (.4) ˆ BG (.4) 式 (.4)~ 式 (.4) を解くことは式 (.6) を解くことと等価である. また 式 (.4)~ 式 (.4) 中に現れる空間的な離散演算子は既に境界条件が考慮されたものであるので 適切に離散化されたフラクショナル ステップ法では中間速度 に対する境界条件は必要ないことがわかる. ^ 次に 行列 B の選択に応じたフラクショナル ステップ法を考える.

123 まず BA - の場合は式 (.5) を直接解くことに等しく分離誤差はゼロであるが A - の計算に時間がかかるのであまり用いられない (Uzw 法 ). KMのフラクショナル ステップ法を含め従来のフラクショナル ステップ法の多くは BI を採用することに対応し その場合 ( AB-I )-.5νL となるので時間 次精度の分離誤差を持つ.KMのフラクショナル ステップ法は一見して時間 次精度のようにも受け取れるが 式 (.7) を分 離して式 (.9)~ 式 (.5) を構成する際に ( ϕ) ( ϕ ) という解析的な互換性が用いられている. しかし 離散式について対応する互換性は一般的に成立しない ( ). GLϕ LGϕ Pero () は時間 次精度および 次精度のフラクショナル ステップ法がそれぞれ B I (.5ν) L および B I (.5ν) L (.5ν) L で構成できることを示した. しかし これらの高次精度のフラクショナル ステップ法はあまり汎用的ではない.

124 4 Dkowz-Dvsky(DD) のフラクショナル ステップ法 (.4) DD DD D G A r r ここで式 (.4) の右辺のベクトルは以下で定義される. b b s r r DD DD D L G ν 非定常問題に対してフラクショナル ステップ法を適用する場合には できれば分離誤差も含めて 次精度以上の時間精度を持つ解法を使用したい.Pero の時間 次精度フラクショナル ステップ法とは別に Dkowz &Dvsky (DD) (6) は式 (.5) をデルタ形式 ( 解の補正量を変数とする形式 ) に書き直した行列方程式を基に より汎用的な時間 次精度のフラクショナル ステップ法が構成できることを示した.

125 DDのフラクショナル ステップ法が時間 次精度となることは式 (.4) の圧力項がデルタ形式であることによる. 速度ベクトルもデルタ形式で記述されているが これは後に式 (.54) で紹介する近似因子化を導入するためであり DDのフラクショナル ステップ法の本質とは関係ない.DDのフラクショナル ステップ法では 式 (.4) の左辺の行列に対してKMのフラクショナル ステップ法と同様な近似を導入する. A D AG r r 式 (.4) から式 (.44) への近似の誤差は AG DD DD (.44) ( ) ( ) ( ) ( G LG O ) ν (.45) となるので DD のフラクショナル ステップ法は時間 次精度の分離誤差を持つ. ここで - O() と評価されている. 5

126 6 (.46) これは次のように分離できる. 近似式 (.44) は次のように LU 分解が可能である. (.47) DD DD I G I DG D A r r DD DD DG D A r r δ δ δ δ I G I (.48)

127 式 (.47) および式 (.48) を展開して整理すると次式を得る. I ν L δ r ˆ δ DGδ D ˆ DD ( ) b ˆ Gδ (.49) (.5) (.5) (.5) δ (.5) 式 (.49)~ 式 (.5) を用いて速度および圧力の時間進行を行うのが DD のフラクショナル ステップ法である. 7

128 DDのフラクショナル ステップ法の計算手順 ) 初期条件として および を与え ) 式 (.49) を解いてδを求め 式 (.5) から^ を計算し ) 式 (.5) を解いてδを求めた後 4) 式 (.5) から を求め 式 (.5) から圧力 を求め 5) および をそれぞれ および として) に戻る. 以上の手順を必要な時間まで繰り返す. DD のフラクショナル ステップ法の計算手順のうち 式 (.49) は速度補正量の離散式 式 (.5) は圧力補正量のポアソン方程式に対応する離散式であり いずれも行列を解く必要がある. 式 (.5) は連続の式の離散式を満足させるために解かれるものであり そのままの形で出来るだけ精度良く解かねばならない. 8

129 式 (.49) に対しては近似因子化 ν I L x ν I L y ν I L z δ r DD (.54) を導入すれば経済的に解くことができる. ここで L x L y L z はそれぞれ x yz 各方向の 階微分に対応する離散演算子であり L と同様に境界条件も考慮されているものとする. 式 (.49) から式 (.54) への演算子の近似誤差は O( ) であるので この近似は式 (.7) が持つ時間精度 ( 次精度 ) を低下させない. 式 (.54) の解 δ は 次元の行列を 回解く事で得られる. ν I ν I ν I L x Ly Lz δ r DD δ δ δ δ (.55) (.55b) (.55) 9

130 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法

131 SIMPLE 法系統の解法 NS 方程式の時間進行に関して MAC 法系統の解法では粘性項および移流項の双方が陽的に時間進行され フラクショナル ステップ法では粘性項が陰的に移流項は陽的に時間進行される.SIMPLE 法系統の解法は粘性項および移流項の双方を陰的に時間進行してさらに大きな時間刻み幅の設定を可能とし収束性および経済性を高める解法である. しかしながら 速度に関して非線形である移流項を完全に陰的に扱うことは各時間ステップの運動方程式の時間積分に繰り返し計算が必要となり経済的でないので SIMPLE 法系統の解法では式 (.) と式 (.) を時間離散化した次式 v v ( v v ) ν v s (.56) v を基に解法が構成される. (.57)

132 (.58) 次元流れでは式 (.56) および式 (.57) は次式となる. (.58b) (.59) x s y x x y v x ν y s y v x v y y v v x v v v ν y v x

133 (.6) さらに式 (.58) および式 (.59) の空間微分を離散化して整理すれば次式を得る. (.6b) (.6) ( ) / / / / / b b b P x b ( ) / / / / / v b b b v P v y b v v / / / / y v v x 式 (.6) 中の総和記号は周辺点 bewns の総和をとることを意味し / に関する式 (.6) の左辺は次式と展開される. ( ) [ ] S S N N W W E E P P b b b P / / / / P E W N S は離散式の係数で 粘性項および移流項の離散化および境界条件に応じて定められる.

134 変数の添字 PEWNSは離散点の相対的な位置を表す. 図.5に / に対するNS 方程式のスタガード格子上の離散点を示す. 速度 / の離散式については P / E / W -/ N / S /- と定義される.SIMPLE 法系統の解法の空間離散化は 通常図中の破線で示されるコントロール ボリュームに関する有限体積法で行われる. 式 (.6) 中の右辺第 項は圧力勾配項の離散式である. 式 (.6) 中のb /は離散式の残りの項をまとめたものである. 速度 v / の離散式も同様に構成され x る ( 図.6). / x / x / x / x / v N x / y / y N y y / v W v P v E y / y / y y y y W P S E y / y y x x v S x x x y / x x x 図.6 v P v / に対する離散点 図.5 P / に対する離散点 4

135 SIMPLE 法 式 (.6) は速度成分 vに関する離散式であるが 右辺に現れる が与えられなければ速度 v は求められない. SIMPLE 法 (7) (Sem-Iml Mehod for Pressre- Lked Eqos) では まず次式を解いて速度の予測値 * v * を求める. ( P ) ( v P ) / / v * / * / b b v b v b * b * b / / b b / v / x y / / (.6) (.6b) しかし * とv * は連続の式の離散式 (.6) を満足しない. そこで圧力および速度の補正を行って式 (.6) を満足させることを考える. δ (.6) * / / δ / v * / v / δv / (.64) (.64b) 5

136 式 (.6) と式 (.64) を式 (.6) に代入した後に式 (.6) との差を計算すると δ および δv と δ との関係式を得る. ( P ) ( v ) P / / δ δv / / b δ b v b δv b b b / / δ δ / (.65) (.65b) ここで 解法が簡単になるように式 (.65) の左辺第 項の総和の項を消去して δ および δv と δ との関係式を簡略化する. x y δ δ / δ δv / / ( P ) x / / δ δ ( v ) P y / / δ δ (.66) (.66b) 6

137 7 (.67) 圧力補正量 δ の方程式は式 (.64)(.66) を連続の式 (.6) に代入して得られる. ここで ( ) ( ) ( ) ( ) ( ) S P v P v P P D y y y x x x / / / / / / / / * δ δ δ δ δ δ δ δ ( ) S y v v x D / / / / * * * * *

138 式 (.67) を解いて得られる圧力補正量 δ を用いて式 (.64) (.66) により速度補正を行うと 連続の式を離散的に満足する速度 および v が得られる. 圧力の補正には式 (.6) ではなく緩和係数 α (<α ) を導入した次式を用いる. α δ (.68) 緩和係数 α の導入は式 (.65) から式 (.66) への簡略化の際に生じる誤差のために必要となる. また 速度場が連続の式を満足するように式 (.67) が定められているので速度の補正に緩和係数は導入しない. 8

139 SIMPLE 法の計算手順 ) 初期条件として v を与え ) 式 (.6) を解いて * と v * を求め ) 式 (.67) を解いて δ を求め 4) 式 (.64) および式 (.66) による速度補正 および式 (.68) による圧力補正を通して v を求め 5) v を v として ) に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. 9

140 なお 式 (.6) のような輸送方程式の離散式を繰り返し法で解く場合 S IMPLE 法系統の解法では減速緩和法を導入して数値計算の安定性を高めることが行われる. 式 (.6) に減速緩和法を導入すれば次式を得る. ( P ) b α / / * / x / b b * b α α / ( P ) re / / (.69) ここで re / は繰り返し計算における つ前のステップでの * の値であり α は減速緩和係数 (<α ) である. 速度 v に対しても また一般的に輸送方程式の離散式を解く場合にも同様の減速緩和法が導入される. SIMPLE 法の問題点は つのパラメータ (α と α ) を流れ場に応じて定めなければならないこと および解の収束性が α の値に大きく依存することである.Pkr (8) は α.5α.8 を推奨しているが その後の研究で α -α とすれば概ね後に述べる SIMPLEC 法と同等の収束性が得られることが示されている (9). 4

141 4 補足 : 行列分解による SIMPLE 法の計算アルゴリズムの導出 (/) 式 (.6) 式 (.6) をまとめて行列表示すれば次式となる. b b D G A ただし G は離散的勾配演算子 D は離散的発散演算子であり これら離散的空間演算子では境界条件も考慮されているものとする. 境界条件から生じる項は NS 方程式および連続の式それぞれに対して b および b に含まれているものとする. ここで δ を代入し 左辺の行列を LU 分解すると以下を得る. b b δ G I G A I G DA D A

142 補足 : 行列分解による SIMPLE 法の計算アルゴリズムの導出 (/) この式は次の 式へと分解できる. A D DA * G δ * b I A G I δ δ G さらに A - に対して次の近似を導入する. [ ( )] A dg A α dg() は行列の対角項を抽出する関数で α は上記対角化近似における補正係数である. * * b 4

143 補足 : 行列分解による SIMPLE 法の計算アルゴリズムの導出 (/) 最後に δ * α δ とおいた後に分解された行列を展開して整理すると次の計算アルゴリズムを得る. A D * b G [ ( )] dg A Gδ D * b [ dg( A) ] Gδ * α δ これは SIMPLE 法の計算アルゴリズムである. 以上の導出から SIMPLE 法の圧力補正係数 α の導入は A - の近似に伴うものであることが理解される. 4

144 SIMPLER 法 SIMPLE 法の収束性を改善する解法として SIMPLE R 法 (SIMPLE Revsed) が提案されている (8).SIM PLE 法の問題点は圧力補正における減速緩和にあるので SIMPLER 法では圧力補正量は連続の式を満足させるためだけに用い 圧力は補正計算をせずに運動方程式から直接求める. この手順を説明するため まず式 (.6) を変形した次式を導入する. v b b b b / / ( P ) ( P ) x / / v v b b v b b / / / v ( P ) ( v P ) y / / / (.7) (.7b) 44

145 式 (.7) の右辺第 項は速度のみで計算される. また解が収束した時点で となることを考慮し 次の中間速度 ^ およびv^ を定義する. ˆ vˆ / / b b v b b b ( P ) / v b v b b v ( P ) / / / (.7) (.7b) この中間速度は既知量である. また 式 (.7) の右辺第 項を^および^v で置き換え この変形に伴う圧力を とすれば次式を得る. / ˆ / (.7) ( ) P x / / v / vˆ / / ( v ) P y / (.7b) 45

146 46 式 (.7) を連続の式の離散式 (.6) に代入すれば圧力 の方程式が得られる. (.7) ここで ( ) ( ) ( ) ( ) ( ) S P v P v P P D y y y x x x / / / / / / / / ˆ ( ) S y v v x D / / / / ˆ ˆ ˆ ˆ ˆ あとは SIMPLE 法と同様に を用いて速度の予測値 * および v * の方程式を解いた後に圧力補正量 δ の方程式を解き速度補正を行う. ただしその際に圧力の補正は行わない.

147 SIMPLER 法の計算手順 ) 初期条件として v を与え ) 式 (.7) から中間速度 ^ およびv^ を計算し ) 式 (.7) を解いて圧力 を求め 4) 式 (.6) を解いて* とv* を求め 5) 式 (.67) を解いてδを求め 6) 式 (.64) および式 (.66) による速度補正を通して速度成分 v を求め 7) v を v として ) に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. SIMPLER 法では圧力方程式を 回 ( 式 (.67) と式 (.7)) 解く必要があるので時間進行 ステップあたりの計算時間はSIMPLE 法よりもかかるが 圧力の緩和係数は使用しないのでその最適値を探す必要がなく 定常計算における収束性は一般的にSIMPLE 法よりも優れている. なお 式 (.6) を繰り返し法で解く過程では式 (.69) と同様の減速緩和法が導入される. 47

148 SIMPLEC 法 SIMPLE 法において圧力の緩和係数を導入しなければならないのは速度補正式の導出の過程における式 (.65) から式 (.66) への大胆な近似に起因する. 速度補正式の導出の過程において SIMPLE 法よりは妥当な近似を用いて速度補正式を導出するのが SIMPLEC 法である (9).SIMPLE C 法の説明のため まず式 (.65) を次の様に書き直す. P b b ( δ δ ) δ δ / b b P / b / x / δ (.74) v P b v b ( δv δv ) δ v δv / b b P / b / y / δ (.74b) 速度補正量は空間に連続的に分布すると考えられるので式 (.74) の左辺第 項は他の項よりも値が小さいことが期待される. 48

149 これより 式 (.74) の左辺第 項を消去して δ および δv と δ との関係式が与えられる. δ δv / / P v P b b b v b / / δ δ x y δ / δ / (.75) (.75b) 49

150 5 圧力補正量 δ の方程式は式 (.64) と式 (.75) を連続の式 (.6) に代入して得られる. (.76) ( ) S b b v P v b b v P v b b P b b P D y y y x x x / / / / / / / / * δ δ δ δ δ δ δ δ 上式を解いて得られる圧力補正量 δ を用いて式 (.64) 式 (.75) により速度補正を行うと連続の式を離散的に満足する速度 および v が得られる. SIMPLE 法とは異なり SIMPLEC 法では圧力の補正に対する緩和計算は必要なく式 (.6) を直接用いる.

151 SIMPLEC 法の計算手順 ) 初期条件として v を与え ) 式 (.6) を解いて * と v* を求め ) 式 (.76) を解いて δ を求め 4) 式 (.64) および式 (.75) による速度補正 および式 (.6) による圧力補正を通して v を求め 5) v を v として ) に戻る. 以上の手順を 非定常計算の場合には必要な時間まで 定常計算の場合には収束条件が満足されるまで繰り返す. SIMPLEC 法では SIMPLER 法と同様に圧力の緩和係数を導入する必要がないことが大きな利点である. また SIMPLEC 法で解かれる圧力補正量の式 (.76) および速度補正量の定義式 (.75) は SI MPLE 法で対応する式の係数が変更されているだけであるので SI MPLEC 法の時間進行 ステップあたりの計算時間は SIMPLE 法とあまり変わらず SIMPLER 法よりは少なくなる. なお 式 (.6) を繰り返し法で解く過程では式 (.69) と同様の減速緩和法が導入される. 5

152 第 章の参考文献 (/) () Hrlow.H. d Welh J.E. Nmerl llo of medeede vsos omressble flow of fld wh free srfe Phys. lds 8 (965.) () Amsde A.A. d Hrlow.H. The SMAC mehod: A merl ehqe for llg omressble fld flows Los Almos sef lborory of he Uv. of Clfor Re. No. LA-47 (97). () Hr C.W. d Cook J.L. Cllg hree-dmesol flows rod srres d over rogh err J. Com. Phys. (97) 4-4. (4) Km J. d Mo P. Alo of frol-se mehod o omressble Nver-Sokes eqos J. Com. Phys. 5.9 (985.) 8-. (5) Pero J.B. A lyss of he frol se mehod J. Com. Phys. 8 (99)

153 第 章の参考文献 (/) (6) Dkowz J.K. d Dvsky A.S. Aroxme forzo s hgh order slg for he ml omressble flow eqos J. Com. Phys. (99) (7) Pker S.V. d Sldg D.B. A llo roedre for he mss d momem rsfer hree-dmesol rbol flows I. J. He Mss Trsfer 5. (97) (8) Pkr S.V. Nmerl he rsfer d fld flow (MGrw-Hll New York) (98) ( 水谷 香月訳 コンピュータによる熱移動と流れの数値解析 ( 森北出版 ) (985.)). (9) V Doorml J.P. d Rhby G.D. Ehemes of he SIMPLE mehod for redg omressble fld flows Nmerl He Trsfer 7 (984) () 日本機械学会 計算力学ハンドブック ( 第 Ⅱ 巻 ) (6) 第. 節 (.-46). 5

154 講義内容 午前 ( 9:-:). 時間進行法 線形多段解法 アダムス バシュフォース公式 アダムス モルトン公式 後退差分公式 ルンゲ クッタ法 陽的ルンゲ クッタ法 陰的ルンゲ クッタ法 時間進行法の安定性解析 線形多段解法の絶対安定領域 ルンゲ クッタ法の絶対安定領域 午後 (:-6:). 非圧縮性流れの計算アルゴリズム 差分格子系と MAC 法系統の解法 レギュラー格子を用いる差分法の問題点 MAC 法とスタガード格子 SMAC 法 4 SOLA 法 (HSMAC 法 ) フラクショナル ステップ法 フラクショナル ステップ法と分離誤差 DD のフラクショナル ステップ法 SIMPLE 法系統の解法 SIMPLE 法 SIMPLER 法 SIMPLEC 法 54